eigen 2.1.4 → 2.1.5
raw patch · 4 files changed
+114/−136 lines, 4 filesdep +binarydep +mtldep ~transformers
Dependencies added: binary, mtl
Dependency ranges changed: transformers
Files
- Data/Eigen/Internal.hsc +28/−23
- Data/Eigen/Matrix.hs +22/−27
- Data/Eigen/SparseMatrix.hs +58/−81
- eigen.cabal +6/−5
Data/Eigen/Internal.hsc view
@@ -1,4 +1,5 @@ {-# OPTIONS_HADDOCK hide #-}+{-# OPTIONS_GHC -fno-warn-orphans #-} {-# LANGUAGE CPP #-} {-# LANGUAGE EmptyDataDecls #-}@@ -22,11 +23,12 @@ #endif import System.IO.Unsafe import Data.Complex-import Data.IORef import Data.Bits+import Data.Binary+import Data.Binary.Get+import Data.Binary.Put import qualified Data.Vector.Storable as VS import qualified Data.ByteString as BS-import qualified Data.ByteString.Lazy as BSL import qualified Data.ByteString.Internal as BSI class (Num a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where@@ -39,6 +41,18 @@ class Cast a b where cast :: a -> b +instance Storable a => Binary (VS.Vector a) where+ put vs = put (BS.length bs) >> putByteString bs where+ (fp,fs) = VS.unsafeToForeignPtr0 vs+ es = sizeOf (VS.head vs)+ bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)+ + get = get >>= getByteString >>= \bs -> let+ (fp,fo,fs) = BSI.toForeignPtr bs+ es = sizeOf (VS.head vs)+ vs = VS.unsafeFromForeignPtr0 (plusForeignPtr fp fo) (fs `div` es)+ in return vs+ -- | Complex number for FFI with the same memory layout as std::complex\<T\> data CComplex a = CComplex !a !a @@ -77,6 +91,9 @@ instance Cast (CComplex CDouble) (Complex Double) where; cast (CComplex x y) = cast x :+ cast y instance Cast (Complex Double) (CComplex CDouble) where; cast (x :+ y) = CComplex (cast x) (cast y) +instance Cast a b => Cast (CTriplet a) (Int, Int, b) where; cast (CTriplet x y z) = (cast x, cast y, cast z)+instance Cast a b => Cast (Int, Int, a) (CTriplet b) where; cast (x,y,z) = CTriplet (cast x) (cast y) (cast z)+ intSize :: Int intSize = sizeOf (undefined :: CInt) @@ -121,9 +138,15 @@ instance Code (CComplex CFloat) where; code _ = 2 instance Code (CComplex CDouble) where; code _ = 3 -magicCode :: Code a => a -> CInt-magicCode x = code x `xor` 0x45696730+newtype MagicCode = MagicCode CInt deriving Eq +instance Binary MagicCode where+ put (MagicCode code) = putWord32be $ fromIntegral code+ get = MagicCode . fromIntegral <$> getWord32be++magicCode :: Code a => a -> MagicCode+magicCode x = MagicCode (code x `xor` 0x45696730)+ #let api1 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name #api1 random, "Ptr b -> CInt -> CInt -> IO CString"@@ -191,6 +214,7 @@ #api2 sparse_setIdentity, "CSparseMatrixPtr a b -> IO CString" #api2 sparse_reserve, "CSparseMatrixPtr a b -> CInt -> IO CString" #api2 sparse_resize, "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"+ #api2 sparse_conservativeResize, "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString" #api2 sparse_compressInplace, "CSparseMatrixPtr a b -> IO CString" #api2 sparse_uncompressInplace, "CSparseMatrixPtr a b -> IO CString"@@ -223,22 +247,3 @@ #api3 sparse_la_logAbsDeterminant, "CSolverPtr a b -> Ptr b -> IO CString" #api3 sparse_la_absDeterminant, "CSolverPtr a b -> Ptr b -> IO CString" #api3 sparse_la_signDeterminant, "CSolverPtr a b -> Ptr b -> IO CString"---openStream :: BSL.ByteString -> IO (IORef BSL.ByteString)-openStream = newIORef--readStream :: IORef BSL.ByteString -> Int -> IO BS.ByteString-readStream ref size = readIORef ref >>= \a ->- let (b,c) = BSL.splitAt (fromIntegral size) a- in if BSL.length b /= fromIntegral size- then fail "readStream: stream exhausted"- else do- writeIORef ref c- return . BS.concat . BSL.toChunks $ b--closeStream :: IORef BSL.ByteString -> IO ()-closeStream ref = BSL.null <$> readIORef ref >>= (`unless` fail "closeStream: stream underrun")--readInt :: IORef BSL.ByteString -> IO CInt-readInt st = decodeInt <$> readStream st intSize
Data/Eigen/Matrix.hs view
@@ -106,8 +106,9 @@ import Prelude hiding (null, sum, all, any, map, filter) import Data.Tuple import Data.Complex hiding (conjugate)+import Data.Binary hiding (encode, decode)+import qualified Data.Binary as B import Foreign.Ptr-import Foreign.ForeignPtr import Foreign.C.Types import Foreign.C.String import Foreign.Storable@@ -125,7 +126,6 @@ import qualified Data.Eigen.Internal as I import qualified Data.Eigen.Matrix.Mutable as M import qualified Data.ByteString.Lazy as BSL-import qualified Data.ByteString.Internal as BSI -- | Matrix to be used in pure computations, uses column major memory layout, features copy-free FFI with C++ <http://eigen.tuxfamily.org Eigen> library. @@ -158,6 +158,26 @@ abs = map abs negate = map negate +-- | Matrix binary serialization+instance I.Elem a b => Binary (Matrix a b) where+ put (Matrix rows cols vals) = do+ put $ I.magicCode (undefined :: b)+ put rows+ put cols+ put vals++ get = do+ get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))+ Matrix <$> get <*> get <*> get++-- | Encode the matrix as a lazy byte string+encode :: I.Elem a b => Matrix a b -> BSL.ByteString+encode = B.encode++-- | Decode matrix from the lazy byte string+decode :: I.Elem a b => BSL.ByteString -> Matrix a b+decode = B.decode+ -- | Empty 0x0 matrix {-# INLINE empty #-} empty :: I.Elem a b => Matrix a b@@ -556,31 +576,6 @@ unsafeWith m@(Matrix rows cols vals) f | not (valid m) = fail "Matrix.unsafeWith: matrix layout is invalid" | otherwise = VS.unsafeWith vals $ \p -> f p (I.cast rows) (I.cast cols)---- | Encode the matrix as a lazy byte string-encode :: I.Elem a b => Matrix a b -> BSL.ByteString-encode m@(Matrix rows cols vals)- | valid m = BSL.fromChunks [- I.encodeInt (I.magicCode $ VS.head vals),- I.encodeInt (I.cast rows),- I.encodeInt (I.cast cols),- let (fp, fs) = VS.unsafeToForeignPtr0 vals in BSI.PS (castForeignPtr fp) 0 (fs * sizeOf (VS.head vals))]- | otherwise = error "Matrix.encode: matrix layout is invalid"---- | Decode matrix from the lazy byte string-decode :: I.Elem a b => BSL.ByteString -> Matrix a b-decode st = Matrix rows cols vals where- (rows, cols, vals) = I.performIO $ do- st <- I.openStream st- code <- I.readInt st- when (code /= I.magicCode (VS.head vals)) $- fail "Matrix.decode: wrong matrix type"- rows <- I.cast <$> I.readInt st- cols <- I.cast <$> I.readInt st- BSI.PS fp fo _ <- I.readStream st (rows * cols * sizeOf (VS.head vals))- I.closeStream st- return (rows, cols, VS.unsafeFromForeignPtr0 (I.plusForeignPtr fp fo) (rows * cols))- {-# INLINE _prop #-} _prop :: I.Elem a b => (Ptr b -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> a
Data/Eigen/SparseMatrix.hs view
@@ -23,6 +23,8 @@ -- * Matrix conversions fromList, toList,+ fromVector,+ toVector, fromDenseList, toDenseList, fromMatrix,@@ -62,6 +64,8 @@ import Prelude hiding (map) import qualified Data.List as L import Data.Complex+import Data.Binary hiding (encode, decode)+import qualified Data.Binary as B import Foreign.C.Types import Foreign.C.String import Foreign.Storable@@ -80,9 +84,7 @@ import qualified Data.Eigen.Internal as I import qualified Data.Vector.Storable as VS import qualified Data.Vector.Storable.Mutable as VSM-import qualified Data.ByteString as BS import qualified Data.ByteString.Lazy as BSL-import qualified Data.ByteString.Internal as BSI {-| A versatible sparse matrix representation. @@ -158,17 +160,30 @@ (+) = add (-) = sub fromInteger x = fromList 1 1 [(0,0,fromInteger x)]- signum = map signum- abs = map abs- negate = map negate+ signum = _map signum+ abs = _map abs+ negate = _map negate --- | Not exposed, For internal use donly-map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b-map f m = fromList (rows m) (cols m) . P.map (\(r,c,v) -> (r,c,f v)) . toList $ m+instance I.Elem a b => Binary (SparseMatrix a b) where+ put m = do+ put $ I.magicCode (undefined :: b)+ put $ rows m+ put $ cols m+ put $ toVector m -mk :: I.Elem a b => Ptr (I.CSparseMatrix a b) -> IO (SparseMatrix a b)-mk p = SparseMatrix <$> FC.newForeignPtr p (I.call $ I.sparse_free p)+ get = do+ get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))+ fromVector <$> get <*> get <*> get +-- | Encode the sparse matrix as a lazy byte string+encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString+encode = B.encode+++-- | Decode sparse matrix from the lazy byte string+decode :: I.Elem a b => BSL.ByteString -> SparseMatrix a b+decode = B.decode+ -- | Stores the coefficient values of the non-zeros. values :: I.Elem a b => SparseMatrix a b -> VS.Vector b values = _getvec I.sparse_values@@ -224,7 +239,7 @@ -- | Extract rectangular block from sparse matrix defined by startRow startCol blockRows blockCols block :: I.Elem a b => Int -> Int -> Int -> Int -> SparseMatrix a b -> SparseMatrix a b-block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) mk+block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) _mk -- | Number of non-zeros elements in the sparse matrix nonZeros :: I.Elem a b => SparseMatrix a b -> Int@@ -232,11 +247,11 @@ -- | The matrix in the compressed format compress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-compress = _unop I.sparse_makeCompressed mk+compress = _unop I.sparse_makeCompressed _mk -- | The matrix in the uncompressed mode uncompress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-uncompress = _unop I.sparse_uncompress mk+uncompress = _unop I.sparse_uncompress _mk -- | Is this in compressed form? compressed :: I.Elem a b => SparseMatrix a b -> Bool@@ -250,51 +265,57 @@ outerSize :: I.Elem a b => SparseMatrix a b -> Int outerSize = _unop I.sparse_outerSize (return . I.cast) --- | Suppresses all nonzeros which are much smaller than reference under the tolerence epsilon+-- | Suppresses all nonzeros which are much smaller than reference under the tolerence @epsilon@ pruned :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b-pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (I.cast r) >> I.sparse_prunedRef p pr pq) mk+pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (I.cast r) >> I.sparse_prunedRef p pr pq) _mk -- | Multiply matrix on a given scalar scale :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b-scale x = _unop (\p pq -> alloca $ \px -> poke px (I.cast x) >> I.sparse_scale p px pq) mk+scale x = _unop (\p pq -> alloca $ \px -> poke px (I.cast x) >> I.sparse_scale p px pq) _mk -- | Transpose of the sparse matrix transpose :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-transpose = _unop I.sparse_transpose mk+transpose = _unop I.sparse_transpose _mk -- | Adjoint of the sparse matrix adjoint :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-adjoint = _unop I.sparse_adjoint mk+adjoint = _unop I.sparse_adjoint _mk -- | Adding two sparse matrices by adding the corresponding entries together. You can use @(+)@ function as well. add :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b-add = _binop I.sparse_add mk+add = _binop I.sparse_add _mk -- | Subtracting two sparse matrices by subtracting the corresponding entries together. You can use @(-)@ function as well. sub :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b-sub = _binop I.sparse_sub mk+sub = _binop I.sparse_sub _mk -- | Matrix multiplication. You can use @(*)@ function as well. mul :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b-mul = _binop I.sparse_mul mk+mul = _binop I.sparse_mul _mk -- | Construct sparse matrix of given size from the list of triplets (row, col, val) fromList :: I.Elem a b => Int -> Int -> [(Int, Int, a)] -> SparseMatrix a b-fromList rows cols tris = I.performIO $ VS.unsafeWith vs $ \p -> alloca $ \pq -> do- I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length vs) pq- peek pq >>= mk- where vs = VS.fromList $ P.map (\(row,col,val) -> I.CTriplet (I.cast row) (I.cast col) (I.cast val)) tris+fromList rows cols = fromVector rows cols . VS.fromList . P.map I.cast +-- | Construct sparse matrix of given size from the storable vector of triplets (row, col, val)+fromVector :: I.Elem a b => Int -> Int -> VS.Vector (I.CTriplet b) -> SparseMatrix a b+fromVector rows cols tris = I.performIO $ VS.unsafeWith tris $ \p -> alloca $ \pq -> do+ I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length tris) pq+ peek pq >>= _mk+ -- | Convert sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included toList :: I.Elem a b => SparseMatrix a b -> [(Int, Int, a)]-toList m@(SparseMatrix fp) = I.performIO $ do+toList = P.map I.cast . VS.toList . toVector++-- | Convert sparse matrix to the storable vector of triplets (row, col, val). Compressed elements will not be included+toVector :: I.Elem a b => SparseMatrix a b -> VS.Vector (I.CTriplet b)+toVector m@(SparseMatrix fp) = I.performIO $ do let size = nonZeros m tris <- VSM.new size withForeignPtr fp $ \p -> VSM.unsafeWith tris $ \q -> I.call $ I.sparse_toList p q (I.cast size)- let f (I.CTriplet row col val) = (I.cast row, I.cast col, I.cast val)- P.map f . VS.toList <$> VS.unsafeFreeze tris+ VS.unsafeFreeze tris -- | Construct sparse matrix of two-dimensional list of values. Matrix dimensions will be detected automatically. Zero values will be compressed. fromDenseList :: (I.Elem a b, Eq a) => [[a]] -> SparseMatrix a b@@ -316,7 +337,7 @@ fromMatrix m1 = I.performIO $ alloca $ \pm0 -> M.unsafeWith m1 $ \vals rows cols -> do I.call $ I.sparse_fromMatrix vals rows cols pm0- peek pm0 >>= mk+ peek pm0 >>= _mk -- | Construct dense matrix from sparse matrix toMatrix :: I.Elem a b => SparseMatrix a b -> M.Matrix a b@@ -327,58 +348,6 @@ I.call $ I.sparse_toMatrix pm1 vals rows cols M.unsafeFreeze m0 --- | Encode the sparse matrix as a lazy byte string-encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString-encode m@(SparseMatrix fp) = I.performIO $ do- let size = nonZeros m- tris <- VSM.new size- withForeignPtr fp $ \p ->- VSM.unsafeWith tris $ \q ->- I.call $ I.sparse_toList p q (I.cast size)- tris <- VS.unsafeFreeze tris- let- tri@(I.CTriplet _ _ val) = VS.head tris-- return $ BSL.fromChunks [- encodeInt (I.magicCode val),- encodeInt (I.cast $ rows m),- encodeInt (I.cast $ cols m),- encodeInt (I.cast $ size),- let (fp, fs) = VS.unsafeToForeignPtr0 tris in BSI.PS (castForeignPtr fp) 0 (fs * sizeOf tri)]-- where- encodeInt :: CInt -> BS.ByteString- encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr----- | Decode sparse matrix from the lazy byte string-decode :: forall a b . I.Elem a b => BSL.ByteString -> SparseMatrix a b-decode st = I.performIO $ do- let- val = undefined :: b- tri = undefined :: I.CTriplet b- triSize = sizeOf tri-- st <- I.openStream st-- code <- I.decodeInt <$> I.readStream st I.intSize- when (code /= I.magicCode val) $- fail "SparseMatrix.decode: wrong matrix type"-- rows <- I.readInt st- cols <- I.readInt st- size <- I.readInt st-- BSI.PS fp fo _ <- I.readStream st (I.cast size * triSize)-- I.closeStream st-- let tris = VS.unsafeFromForeignPtr0 (I.plusForeignPtr fp fo) (I.cast size)-- VS.unsafeWith tris $ \p -> alloca $ \pq -> do- I.call $ I.sparse_fromList rows cols p size pq- peek pq >>= mk- -- | Yield an immutable copy of the mutable matrix freeze :: I.Elem a b => SMM.IOSparseMatrix a b -> IO (SparseMatrix a b) freeze (SMM.IOSparseMatrix fp) = SparseMatrix <$> _clone fp@@ -426,3 +395,11 @@ I.call $ I.sparse_clone p pq q <- peek pq FC.newForeignPtr q $ I.call $ I.sparse_free q++_map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b+_map f m = fromVector (rows m) (cols m) . VS.map g . toVector $ m where+ g (I.CTriplet r c v) = I.CTriplet r c $ I.cast $ f $ I.cast v++_mk :: I.Elem a b => Ptr (I.CSparseMatrix a b) -> IO (SparseMatrix a b)+_mk p = SparseMatrix <$> FC.newForeignPtr p (I.call $ I.sparse_free p)+
eigen.cabal view
@@ -1,5 +1,5 @@ name: eigen-version: 2.1.4+version: 2.1.5 homepage: https://github.com/osidorkin/haskell-eigen synopsis: Eigen C++ library (linear algebra: matrices, sparse matrices, vectors, numerical solvers). description: This module provides Haskell binding for <http://eigen.tuxfamily.org/ Eigen C++ library>.@@ -1308,6 +1308,7 @@ build-depends: base >= 3 && < 5, vector >= 0.5 && < 0.12, primitive >= 0.1 && < 0.7,+ binary, bytestring, transformers >= 0.3 @@ -1322,20 +1323,20 @@ Test-Suite test-solve type: exitcode-stdio-1.0 main-is: test/solve.hs- build-depends: base, primitive, vector, eigen+ build-depends: base, primitive, vector, bytestring, transformers, binary, eigen Test-Suite test-solve-sparse type: exitcode-stdio-1.0 main-is: test/solve-sparse.hs- build-depends: base, primitive, vector, eigen+ build-depends: base, primitive, vector, bytestring, transformers, mtl, binary, eigen Test-Suite test-rank type: exitcode-stdio-1.0 main-is: test/rank.hs- build-depends: base, primitive, vector, eigen+ build-depends: base, primitive, vector, bytestring, transformers, binary, eigen Test-Suite test-regression type: exitcode-stdio-1.0 main-is: test/regression.hs- build-depends: base, primitive, vector, eigen+ build-depends: base, primitive, vector, bytestring, transformers, binary, eigen