hmatrix 0.10.0.2 → 0.11.0.0

raw patch · 12 files changed

+409/−31 lines, 12 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• CHANGES +13/−0
• THANKS +3/−0
• examples/bool.hs +51/−0
• hmatrix.cabal +15/−2
• lib/Data/Packed/Internal/Matrix.hs +33/−1
• lib/Data/Packed/Internal/Vector.hs +35/−0
• lib/Data/Packed/Matrix.hs +10/−18
• lib/Numeric/ContainerBoot.hs +122/−0
• lib/Numeric/LinearAlgebra/Algorithms.hs +20/−5
• lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +42/−0
• lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h +6/−0
• lib/Numeric/LinearAlgebra/Tests.hs +59/−5

@@ -1,3 +1,16 @@+0.11.0.0+========++- flag -fvector default = True++- invlndet (inverse and log of determinant)++- step, cond++- find++- assoc, accum+ 0.10.0.0 ======== 

@@ -81,3 +81,6 @@  - Alexey Khudyakov cleaned up PRAGMAS and fixed some hlint suggestions. +- Torsten Kemps-Benedix reported an installation problem in OS/X.++

@@ -0,0 +1,51 @@+-- vectorized boolean operations defined in terms of step or cond++import Numeric.LinearAlgebra++infix  4  .==., ./=., .<., .<=., .>=., .>.+infixr 3  .&&.+infixr 2  .||.++a .<.  b = step (b-a)+a .<=. b = cond a b 1 1 0+a .==. b = cond a b 0 1 0+a ./=. b = cond a b 1 0 1+a .>=. b = cond a b 0 1 1+a .>.  b = step (a-b)++a .&&. b  = step (a*b)+a .||. b  = step (a+b)+no a      = 1-a+xor a b   = a ./=. b+equiv a b = a .==. b+imp a b   = no a .||. b++taut x = minElement x == 1++-- examples++clip a b x = cond y b y y b where y = cond x a a x x++disp = putStr . dispf 3++eye n = ident n :: Matrix Double+row = asRow . fromList    :: [Double] -> Matrix Double+col = asColumn . fromList :: [Double] -> Matrix Double++m = (3><4) [1..] :: Matrix Double++p = row [0,0,1,1]+q = row [0,1,0,1]++main = do+    print $ find (>6) m+    disp $ assoc (6,8) 7 $ zip (find (/=0) (eye 5)) [10..]+    disp $ accum (eye 5) (+) [((0,2),3), ((3,1),7), ((1,1),1)]+    disp $ m .>=. 10  .||.  m .<. 4+    (disp . fromColumns . map flatten) [p, q, p.&&.q, p .||.q, p `xor` q, p `equiv` q, p `imp` q]+    print $ taut $ (p `imp` q ) `equiv` (no q `imp` no p)+    print $ taut $ (xor p q) `equiv` (p .&&. no q .||. no p .&&. q)+    disp $ clip 3 8 m+    disp $ col [1..7] .<=. row [1..5]+    disp $ cond (col [1..3]) (row [1..4]) m 50 (3*m)+

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.10.0.2+Version:            0.11.0.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -54,6 +54,7 @@                     examples/Real.hs                     examples/vector.hs                     examples/monadic.hs+                    examples/bool.hs  extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,                     lib/Numeric/LinearAlgebra/LAPACK/clapack.h@@ -76,12 +77,16 @@  flag vector     description:    Use Data.Vector.Storable type from "vector" package.-    default:        False+    default:        True  flag binary     description:    Define Binary instances     default:        True +flag finit+    description:    Force FPU initialization in foreing calls+    default:        False+ library      Build-Depends:      base >= 4 && < 5,@@ -159,6 +164,13 @@     if impl(ghc < 6.10.2)         cpp-options: -DFINIT +    if flag(finit)+        cpp-options: -DFINIT++    if impl(ghc == 7.0.1)+        cpp-options: -DNONORMVTEST+        cpp-options: -DFINIT+     if flag(mkl)       if arch(x86_64)         extra-libraries:   gsl mkl_lapack mkl_intel_lp64 mkl_sequential mkl_core@@ -169,6 +181,7 @@         extra-lib-dirs: /opt/local/lib/         include-dirs: /opt/local/include         extra-libraries: gsl+        cc-options: -arch i386         frameworks: Accelerate      if os(windows)

@@ -31,7 +31,8 @@     liftMatrix, liftMatrix2,     (@@>),     saveMatrix,-    singleton+    singleton,+    size, shSize, conformVs, conformMs, conformVTo, conformMTo ) where  import Data.Packed.Internal.Common@@ -441,3 +442,34 @@     free charfmt  foreign import ccall "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM++----------------------------------------------------------------------++conformMs ms = map (conformMTo (r,c)) ms+  where+    r = maximum (map rows ms)+    c = maximum (map cols ms)++conformVs vs = map (conformVTo n) vs+  where+    n = maximum (map dim vs)++conformMTo (r,c) m+    | size m == (r,c) = m+    | size m == (1,1) = reshape c (constantD (m@@>(0,0)) (r*c))+    | size m == (r,1) = repCols c m+    | size m == (1,c) = repRows r m+    | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"++conformVTo n v+    | dim v == n = v+    | dim v == 1 = constantD (v@>0) n+    | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n++repRows n x = fromRows (replicate n (flatten x))+repCols n x = fromColumns (replicate n (flatten x))++size m = (rows m, cols m)++shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"+

@@ -22,6 +22,7 @@     foldVector, foldVectorG, foldLoop, foldVectorWithIndex,     createVector, vec,     asComplex, asReal, float2DoubleV, double2FloatV,+    stepF, stepD, condF, condD,     fwriteVector, freadVector, fprintfVector, fscanfVector,     cloneVector,     unsafeToForeignPtr,@@ -291,6 +292,40 @@  foreign import ccall "float2double" c_float2double:: TFV foreign import ccall "double2float" c_double2float:: TVF++---------------------------------------------------------------++stepF :: Vector Float -> Vector Float+stepF v = unsafePerformIO $ do+    r <- createVector (dim v)+    app2 c_stepF vec v vec r "stepF"+    return r++stepD :: Vector Double -> Vector Double+stepD v = unsafePerformIO $ do+    r <- createVector (dim v)+    app2 c_stepD vec v vec r "stepD"+    return r++foreign import ccall "stepF" c_stepF :: TFF+foreign import ccall "stepD" c_stepD :: TVV++---------------------------------------------------------------++condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float+condF x y l e g = unsafePerformIO $ do+    r <- createVector (dim x)+    app6 c_condF vec x vec y vec l vec e vec g vec r "condF"+    return r++condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double+condD x y l e g = unsafePerformIO $ do+    r <- createVector (dim x)+    app6 c_condD vec x vec y vec l vec e vec g vec r "condD"+    return r++foreign import ccall "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF+foreign import ccall "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV  ---------------------------------------------------------------- 

@@ -302,29 +302,21 @@ -- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix. liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t liftMatrix2Auto f m1 m2-    | compat' m1 m2 = lM f m1 m2--    | r1 == 1  && c2 == 1 = lM f (repRows r2 m1) (repCols c1 m2)-    | c1 == 1  && r2 == 1 = lM f (repCols c2 m1) (repRows r1 m2)-    -    | r1 == r2 && c2 == 1 = lM f m1              (repCols c1 m2)-    | r1 == r2 && c1 == 1 = lM f (repCols c2 m1) m2--    | c1 == c2 && r2 == 1 = lM f m1              (repRows r1 m2)-    | c1 == c2 && r1 == 1 = lM f (repRows r2 m1) m2--    | otherwise    = error $ "nonconformable matrices in liftMatrix2Auto: "-                             ++ show (size m1) ++ ", " ++ show (size m2)+    | compat' m1 m2 = lM f m1  m2+    | ok            = lM f m1' m2'+    | otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ shSize m1 ++ ", " ++ shSize m2   where     (r1,c1) = size m1     (r2,c2) = size m2--size m = (rows m, cols m)+    r = max r1 r2+    c = max c1 c2+    r0 = min r1 r2+    c0 = min c1 c2+    ok = r0 == 1 || r1 == r2 && c0 == 1 || c1 == c2+    m1' = conformMTo (r,c) m1+    m2' = conformMTo (r,c) m2  lM f m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))--repRows n x = fromRows (replicate n (flatten x))-repCols n x = fromColumns (replicate n (flatten x))  compat' :: Matrix a -> Matrix b -> Bool compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2

@@ -45,6 +45,7 @@ ) where  import Data.Packed+import Data.Packed.ST as ST import Numeric.Conversion import Data.Packed.Internal import Numeric.GSL.Vector@@ -121,6 +122,63 @@     -- | the product of elements (faster than using @fold@)     prodElements :: c e -> e +    -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@+    --+    -- @> step $ linspace 5 (-1,1::Double)+    -- 5 |> [0.0,0.0,0.0,1.0,1.0]@+    +    step :: RealElement e => c e -> c e++    -- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@.+    --+    -- Arguments with any dimension = 1 are automatically expanded: +    --+    -- @> cond ((1>\<4)[1..]) ((3>\<1)[1..]) 0 100 ((3>\<4)[1..]) :: Matrix Double+    -- (3><4)+    -- [ 100.0,   2.0,   3.0,  4.0+    -- ,   0.0, 100.0,   7.0,  8.0+    -- ,   0.0,   0.0, 100.0, 12.0 ]@+    +    cond :: RealElement e +         => c e -- ^ a+         -> c e -- ^ b+         -> c e -- ^ l +         -> c e -- ^ e+         -> c e -- ^ g+         -> c e -- ^ result++    -- | Find index of elements which satisfy a predicate+    --+    -- @> find (>0) (ident 3 :: Matrix Double)+    -- [(0,0),(1,1),(2,2)]@++    find :: (e -> Bool) -> c e -> [IndexOf c]++    -- | Create a structure from an association list+    --+    -- @> assoc 5 0 [(2,7),(1,3)] :: Vector Double+    -- 5 |> [0.0,3.0,7.0,0.0,0.0]@+    +    assoc :: IndexOf c        -- ^ size+          -> e                -- ^ default value+          -> [(IndexOf c, e)] -- ^ association list+          -> c e              -- ^ result++    -- | Modify a structure using an update function+    --+    -- @> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double+    -- (5><5)+    --  [ 1.0, 0.0, 0.0, 3.0, 0.0+    --  , 0.0, 6.0, 0.0, 0.0, 0.0+    --  , 0.0, 0.0, 1.0, 0.0, 0.0+    --  , 0.0, 0.0, 0.0, 1.0, 0.0+    --  , 0.0, 0.0, 0.0, 0.0, 1.0 ]@+    +    accum :: c e              -- ^ initial structure+          -> (e -> e -> e)    -- ^ update function+          -> [(IndexOf c, e)] -- ^ association list+          -> c e              -- ^ result+ --------------------------------------------------------------------------  instance Container Vector Float where@@ -145,6 +203,11 @@     maxElement  = toScalarF Max     sumElements  = sumF     prodElements = prodF+    step = stepF+    find = findV+    assoc = assocV+    accum = accumV+    cond = condV condF  instance Container Vector Double where     scale = vectorMapValR Scale@@ -168,6 +231,11 @@     maxElement  = toScalarR Max     sumElements  = sumR     prodElements = prodR+    step = stepD+    find = findV+    assoc = assocV+    accum = accumV+    cond = condV condD  instance Container Vector (Complex Double) where     scale = vectorMapValC Scale@@ -191,6 +259,11 @@     maxElement  = ap (@>) maxIndex     sumElements  = sumC     prodElements = prodC+    step = undefined -- cannot match+    find = findV+    assoc = assocV+    accum = accumV+    cond = undefined -- cannot match  instance Container Vector (Complex Float) where     scale = vectorMapValQ Scale@@ -214,6 +287,11 @@     maxElement  = ap (@>) maxIndex     sumElements  = sumQ     prodElements = prodQ+    step = undefined -- cannot match+    find = findV+    assoc = assocV+    accum = accumV+    cond = undefined -- cannot match  --------------------------------------------------------------- @@ -243,6 +321,11 @@     maxElement = ap (@@>) maxIndex     sumElements = sumElements . flatten     prodElements = prodElements . flatten+    step = liftMatrix step+    find = findM+    assoc = assocM+    accum = accumM+    cond = condM  ---------------------------------------------------- @@ -580,3 +663,42 @@ -- | creates the identity matrix of given dimension ident :: (Num a, Element a) => Int -> Matrix a ident n = diag (constantD 1 n)++--------------------------------------------------------++findV p x = foldVectorWithIndex g [] x where+    g k z l = if p z then k:l else l++findM p x = map ((`divMod` cols x)) $ findV p (flatten x)++assocV n z xs = ST.runSTVector $ do+        v <- ST.newVector z n+        mapM_ (\(k,x) -> ST.writeVector v k x) xs+        return v++assocM (r,c) z xs = ST.runSTMatrix $ do+        m <- ST.newMatrix z r c+        mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs+        return m++accumV v0 f xs = ST.runSTVector $ do+        v <- ST.thawVector v0+        mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs+        return v++accumM m0 f xs = ST.runSTMatrix $ do+        m <- ST.thawMatrix m0+        mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs+        return m++----------------------------------------------------------------------++condM a b l e t = reshape (cols a'') $ cond a' b' l' e' t'+  where+    args@(a'':_) = conformMs [a,b,l,e,t]+    [a', b', l', e', t'] = map flatten args++condV f a b l e t = f a' b' l' e' t'+  where+    [a', b', l', e', t'] = conformVs [a,b,l,e,t]+

@@ -31,7 +31,8 @@     linearSolveLS,     linearSolveSVD,     inv, pinv,-    det, rank, rcond,+    det, invlndet,+    rank, rcond, -- * Matrix factorizations -- ** Singular value decomposition     svd,@@ -341,12 +342,25 @@        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"  +-- | Joint computation of inverse and logarithm of determinant of a square matrix.+invlndet :: (Floating t, Field t)+         => Matrix t+         -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) +invlndet m | square m = (im,(ladm,sdm))+           | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"+  where+    lp@(lup,perm) = luPacked m+    s = signlp (rows m) perm+    dg = toList $ takeDiag $ lup+    ladm = sum $ map (log.abs) dg+    sdm = s* product (map signum dg)+    im = luSolve lp (ident (rows m))  --- | Determinant of a square matrix.+-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'. det :: Field t => Matrix t -> t det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)-      | otherwise = error "det of nonsquare matrix"+      | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"     where (lup,perm) = luPacked m           s = signlp (rows m) perm @@ -357,10 +371,10 @@ lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) lu = luFact . luPacked --- | Inverse of a square matrix.+-- | Inverse of a square matrix. See also 'invlndet'. inv :: Field t => Matrix t -> Matrix t inv m | square m = m `linearSolve` ident (rows m)-      | otherwise = error "inv of nonsquare matrix"+      | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"  -- | Pseudoinverse of a general matrix. pinv :: Field t => Matrix t -> Matrix t@@ -697,3 +711,4 @@ relativeError x y = dig (norm (x `sub` y) / norm x)     where norm = pnorm Infinity           dig r = round $ -logBase 10 (realToFrac r :: Double)+

@@ -1247,3 +1247,45 @@     OK } +//////////////////// step /////////////////////////++int stepF(FVEC(x),FVEC(y)) {+    DEBUGMSG("stepF")+    int k;+    for(k=0;k<xn;k++) {+        yp[k]=xp[k]>0;+    }+    OK+}++int stepD(DVEC(x),DVEC(y)) {+    DEBUGMSG("stepD")+    int k;+    for(k=0;k<xn;k++) {+        yp[k]=xp[k]>0;+    }+    OK+}++//////////////////// cond /////////////////////////++int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) {+    REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);+    DEBUGMSG("condF")+    int k;+    for(k=0;k<xn;k++) {+        rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);+    }+    OK+}++int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) {+    REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE);+    DEBUGMSG("condD")+    int k;+    for(k=0;k<xn;k++) {+        rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]);+    }+    OK+}+

@@ -87,6 +87,12 @@ int conjugateQ(KQVEC(x),QVEC(t)); int conjugateC(KCVEC(x),CVEC(t)); +int stepF(FVEC(x),FVEC(y));+int stepD(DVEC(x),DVEC(y));++int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r));+int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r));+ int svd_l_R(KDMAT(x),DMAT(u),DVEC(s),DMAT(v)); int svd_l_Rdd(KDMAT(x),DMAT(u),DVEC(s),DMAT(v)); int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));

@@ -67,6 +67,14 @@             , 2, 8, i             ] +detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]+  where+    m = complex (feye 6)+    inv1 = inv m+    det1 = det m+    (inv2,(lda,sa)) = invlndet m+    det2 = sa * exp lda+ --------------------------------------------------------------------  polyEval cs x = foldr (\c ac->ac*x+c) 0 cs@@ -224,9 +232,10 @@ normsVTest = TestList [     utest "normv2CD" $ norm2PropC v   , utest "normv2CF" $ norm2PropC (single v)+#ifndef NONORMVTEST   , utest "normv2D"  $ norm2PropR x   , utest "normv2F"  $ norm2PropR (single x)-+#endif   , utest "normv1CD" $ norm1 v          == 8   , utest "normv1CF" $ norm1 (single v) == 8   , utest "normv1D"  $ norm1 x          == 6@@ -239,7 +248,9 @@   ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)          x = fromList [1,2,-3] :: Vector Double+#ifndef NONORMVTEST          norm2PropR a = norm2 a =~= sqrt (dot a a)+#endif          norm2PropC a = norm2 a =~= realPart (sqrt (dot a (conj a)))          a =~= b = fromList [a] |~| fromList [b] @@ -334,8 +345,8 @@  -- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs --successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool-successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)-   where step e = do+successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0)+   where stp e  = do                   ep <- lift_maybe $ state_get                   if t e ep                      then lift_maybe $ state_put e@@ -343,8 +354,8 @@  -- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input --successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b-successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)-   where step e = do+successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0)+   where stp  e = do                   ep <- state_get                   state_put e                   return $ f ep e@@ -357,7 +368,45 @@  --------------------------------------------------------------------- +findAssocTest = utest "findAssoc" ok+  where+    ok = m1 == m2+    m1 = assoc (6,6) 7 $ zip (find (>0) (ident 5 :: Matrix Float)) [10 ..] :: Matrix Double+    m2 = diagRect 7 (fromList[10..14]) 6 6 +---------------------------------------------------------------------++condTest = utest "cond" ok+  where+    ok = step v * v == cond v 0 0 0 v+    v = fromList [-7 .. 7 ] :: Vector Float++---------------------------------------------------------------------++conformTest = utest "conform" ok+  where+    ok = 1 + row [1,2,3] + col [10,20,30,40] + (4><3) [1..]+         == (4><3) [13,15,17+                   ,26,28,30+                   ,39,41,43+                   ,52,54,56]+    row = asRow . fromList+    col = asColumn . fromList :: [Double] -> Matrix Double++---------------------------------------------------------------------++accumTest = utest "accum" ok+  where+    x = ident 3 :: Matrix Double+    ok = accum x (+) [((1,2),7), ((2,2),3)]+         == (3><3) [1,0,0+                   ,0,1,7+                   ,0,0,4]+         &&+         toList (flatten x) == [1,0,0,0,1,0,0,0,1] ++---------------------------------------------------------------------+ -- | All tests must pass with a maximum dimension of about 20 --  (some tests may fail with bigger sizes due to precision loss). runTests :: Int  -- ^ maximum dimension@@ -493,6 +542,7 @@     _ <- runTestTT $ TestList         [ utest "1E5 rots" rotTest         , utest "det1" detTest1+        , utest "invlndet" detTest2         , utest "expm1" (expmTest1)         , utest "expm2" (expmTest2)         , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)@@ -524,6 +574,10 @@         , sumprodTest         , chainTest         , succTest+        , findAssocTest+        , condTest+        , conformTest+        , accumTest         ]     return ()