eigen 1.2.2 → 1.2.3

raw patch · 4 files changed

+95/−62 lines, 4 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• Data/Eigen/Internal.hs +6/−0
• Data/Eigen/LA.hs +81/−52
• Data/Eigen/Matrix.hs +5/−8
• eigen.cabal +3/−2

@@ -54,3 +54,9 @@ foreign import ccall "eigen-proxy.h eigen_hypotNorm"   c_hypotNorm   :: Ptr CDouble -> CInt -> CInt -> IO CDouble foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: Ptr CDouble -> CInt -> CInt -> IO CDouble +foreign import ccall "eigen-proxy.h eigen_rank"         c_rank       :: CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_image"        c_image      :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_kernel"       c_kernel     :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h free"               c_free       :: Ptr a -> IO ()+foreign import ccall "eigen-proxy.h eigen_solve"        c_solve       :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString

@@ -48,25 +48,63 @@  Checking if a solution really exists: Only you know what error margin you want to allow for a solution to be considered valid. -You can compute relative error using @norm (ax - b) / norm b@ formula or use 'relativeError' function which provides the same calculation implemented slightly more efficient.+You can compute relative error using @'norm' (ax - b) / 'norm' b@ formula or use 'relativeError' function which provides the same calculation implemented slightly more efficient.  -}  module Data.Eigen.LA (     -- * Basic linear solving     Decomposition(..),+    solve,+    relativeError,+    -- * Rank-revealing decompositions+    {- | +Certain decompositions are rank-revealing, i.e. are able to compute the 'rank' of a matrix. These are typically also the decompositions that behave best in the face of a non-full-rank matrix (which in the 'square' case means a singular matrix).++@+import Data.Eigen.Matrix+import Data.Eigen.LA++main = do+    let a = fromList [[1,2,5],[2,1,4],[3,0,3]]+    putStrLn "Here is the matrix A:" >> print a+    putStrLn "The rank of A is:" >> print (rank FullPivLU a)+    putStrLn "Here is a matrix whose columns form a basis of the null-space of A:" >> print (kernel FullPivLU a)+    putStrLn "Here is a matrix whose columns form a basis of the column-space of A:" >> print (image FullPivLU a)+@++produces the following output++@+Here is the matrix A:+Matrix 3x3+1.0 2.0 5.0+2.0 1.0 4.0+3.0 0.0 3.0++The rank of A is:+2+Here is a matrix whose columns form a basis of the null-space of A:+Matrix 3x1+0.5000000000000001+1.0+-0.5++Here is a matrix whose columns form a basis of the column-space of A:+Matrix 3x2+5.0 1.0+4.0 2.0+3.0 3.0+@+    -}     rank,     kernel,     image,-    solve,-    relativeError,     -- * Multiple linear regression+    {- | A linear regression model that contains more than one predictor variable. -}     linearRegression ) where -import Foreign.Ptr-import Foreign.C.Types-import Foreign.C.String import Foreign.Storable import Foreign.Marshal.Alloc import qualified Foreign.Concurrent as FC@@ -76,15 +114,6 @@ import qualified Data.Eigen.Matrix.Mutable as M import qualified Data.Vector.Storable as VS -foreign import ccall "eigen-proxy.h eigen_solve"    c_solve  :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_rank"     c_rank   :: CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_image"    c_image  :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_kernel"   c_kernel :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h free"           c_free   :: Ptr a -> IO ()--- {- | @ Decomposition           Requirements on the matrix          Speed   Accuracy  Rank  Kernel  Image @@ -97,9 +126,8 @@ LLT                     Positive definite                   +++     +         -     -       - LDLT                    Positive or negative semidefinite   +++     ++        -     -       - JacobiSVD               None                                -       +++       +     -       ---The best way to do least squares solving for square matrices is with a SVD decomposition (JacobiSVD) @+The best way to do least squares solving for square matrices is with a SVD decomposition ('JacobiSVD') -}  data Decomposition@@ -120,12 +148,38 @@     -- | Two-sided Jacobi SVD decomposition of a rectangular matrix.     | JacobiSVD deriving (Show, Enum) ++-- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@+solve :: Decomposition -> Matrix -> Matrix -> Matrix+solve d a b = I.performIO $ do+    x <- M.new (m_cols a) 1+    M.unsafeWith x $ \x_vals x_rows x_cols ->+        unsafeWith a $ \a_vals a_rows a_cols ->+            unsafeWith b $ \b_vals b_rows b_cols ->+                I.call $ I.c_solve (I.cast $ fromEnum d)+                    x_vals x_rows x_cols+                    a_vals a_rows a_cols+                    b_vals b_rows b_cols+    unsafeFreeze x++-- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm+relativeError :: Matrix -> Matrix -> Matrix -> Double+relativeError x a b = I.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+                    I.call $ I.c_relativeError pe+                        x_vals x_rows x_cols+                        a_vals a_rows a_cols+                        b_vals b_rows b_cols+                    I.cast <$> peek pe++-- | The rank of the matrix rank :: Decomposition -> Matrix -> Int-rank d m = I.performIO $-    unsafeWith m $ \vals rows cols ->-        alloca $ \pr -> do-            I.call $ c_rank (I.cast $ fromEnum d) pr vals rows cols-            I.cast <$> peek pr+rank d m = I.performIO $ alloca $ \pr -> do+    I.call $ unsafeWith m $ I.c_rank (I.cast $ fromEnum d) pr+    I.cast <$> peek pr  -- | Return matrix whose columns form a basis of the null-space of @A@ kernel :: Decomposition -> Matrix -> Matrix@@ -134,13 +188,13 @@     alloca $ \prows ->     alloca $ \pcols ->         unsafeWith m1 $ \vals1 rows1 cols1 -> do-            I.call $ c_kernel (I.cast $ fromEnum d)+            I.call $ I.c_kernel (I.cast $ fromEnum d)                 pvals prows pcols                 vals1 rows1 cols1             vals <- peek pvals             rows <- I.cast <$> peek prows             cols <- I.cast <$> peek pcols-            fp <- FC.newForeignPtr vals $ c_free vals+            fp <- FC.newForeignPtr vals $ I.c_free vals             return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols  @@ -151,41 +205,16 @@     alloca $ \prows ->     alloca $ \pcols ->         unsafeWith m1 $ \vals1 rows1 cols1 -> do-            I.call $ c_image (I.cast $ fromEnum d)+            I.call $ I.c_image (I.cast $ fromEnum d)                 pvals prows pcols                 vals1 rows1 cols1             vals <- peek pvals             rows <- I.cast <$> peek prows             cols <- I.cast <$> peek pcols-            fp <- FC.newForeignPtr vals $ c_free vals+            fp <- FC.newForeignPtr vals $ I.c_free vals             return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols --- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@-solve :: Decomposition -> Matrix -> Matrix -> Matrix-solve d a b = I.performIO $ do-    x <- M.new (m_cols a) 1-    M.unsafeWith x $ \x_vals x_rows x_cols ->-        unsafeWith a $ \a_vals a_rows a_cols ->-            unsafeWith b $ \b_vals b_rows b_cols ->-                I.call $ c_solve (I.cast $ fromEnum d)-                    x_vals x_rows x_cols-                    a_vals a_rows a_cols-                    b_vals b_rows b_cols-    unsafeFreeze x --- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm-relativeError :: Matrix -> Matrix -> Matrix -> Double-relativeError x a b = I.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-                    I.call $ c_relativeError pe-                        x_vals x_rows x_cols-                        a_vals a_rows a_cols-                        b_vals b_rows b_cols-                    I.cast <$> peek pe- {- | [(coeffs, error) = linearRegression points] computes multiple linear regression @y = a1 x1 + a2 x2 + ... + an xn + b@ using 'ColPivHouseholderQR' decomposition @@ -203,7 +232,7 @@     [-4.01, 2.41, 6.01],     [-3.86, 2.09, 5.55],     [-4.10, 2.58, 6.32]]- @+@   produces the following output 

@@ -253,11 +253,11 @@ trace :: Matrix -> Double trace = _prop I.c_trace --- | Returns true if all of the coefficients in a given matrix evaluate to true+-- | Applied to a predicate and a matrix, all determines if all elements of the matrix satisfies the predicate all :: (Double -> Bool) -> Matrix -> Bool all f = VS.all (f . I.cast) . m_vals --- | Returns true if at least one of the coefficients in a given matrix evaluates to true+-- | Applied to a predicate and a matrix, any determines if any element of the matrix satisfies the predicate any :: (Double -> Bool) -> Matrix -> Bool any f = VS.any (f . I.cast) . m_vals @@ -287,22 +287,19 @@     | square m = _prop I.c_determinant m     | otherwise = error "Matrix.determinant: non-square matrix" --- | Adding two matrices by adding the corresponding entries together+-- | Adding two matrices by adding the corresponding entries together. You can use @(+)@ function as well. 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" --- | Return a + b+-- | Subtracting two matrices by subtracting the corresponding entries together. You can use @(-)@ function as well. 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" -{- | Matrix multiplication--<<http://upload.wikimedia.org/math/7/3/f/73fc7ef1bf6f6822115c41cff58535e1.png>>--}+-- | Matrix multiplication. You can use @(*)@ function as well. mul :: Matrix -> Matrix -> Matrix mul m1 m2     | m_cols m1 == m_rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.c_mul m1 m2

@@ -1,5 +1,5 @@ name:           eigen-version:        1.2.2+version:        1.2.3 homepage:       https://github.com/osidorkin/haskell-eigen synopsis:       Eigen C++ library (linear algebra: matrices, vectors, numerical solvers). description:    This module provides Haskell binding for Eigen C++ library.@@ -365,7 +365,8 @@                         Data.Eigen.Matrix                         Data.Eigen.Matrix.Mutable                         Data.Eigen.Parallel-                        Data.Eigen.Internal+                        +    other-modules:      Data.Eigen.Internal      ghc-options:        -Wall -fno-warn-name-shadowing     build-depends:      base >= 3 && < 5,