summaryrefslogtreecommitdiff
path: root/packages/base/src/Data/Packed
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Data/Packed')
-rw-r--r--packages/base/src/Data/Packed/Development.hs32
-rw-r--r--packages/base/src/Data/Packed/Foreign.hs100
-rw-r--r--packages/base/src/Data/Packed/IO.hs167
-rw-r--r--packages/base/src/Data/Packed/Internal.hs24
-rw-r--r--packages/base/src/Data/Packed/Internal/Common.hs160
-rw-r--r--packages/base/src/Data/Packed/Internal/Matrix.hs423
-rw-r--r--packages/base/src/Data/Packed/Internal/Numeric.hs720
-rw-r--r--packages/base/src/Data/Packed/Internal/Signatures.hs70
-rw-r--r--packages/base/src/Data/Packed/Internal/Vector.hs471
-rw-r--r--packages/base/src/Data/Packed/Matrix.hs494
-rw-r--r--packages/base/src/Data/Packed/Numeric.hs299
-rw-r--r--packages/base/src/Data/Packed/ST.hs178
-rw-r--r--packages/base/src/Data/Packed/Vector.hs125
13 files changed, 0 insertions, 3263 deletions
diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs
deleted file mode 100644
index 72eb16b..0000000
--- a/packages/base/src/Data/Packed/Development.hs
+++ /dev/null
@@ -1,32 +0,0 @@
1
2-----------------------------------------------------------------------------
3-- |
4-- Module : Data.Packed.Development
5-- Copyright : (c) Alberto Ruiz 2009
6-- License : BSD3
7-- Maintainer : Alberto Ruiz
8-- Stability : provisional
9-- Portability : portable
10--
11-- The library can be easily extended with additional foreign functions
12-- using the tools in this module. Illustrative usage examples can be found
13-- in the @examples\/devel@ folder included in the package.
14--
15-----------------------------------------------------------------------------
16{-# OPTIONS_HADDOCK hide #-}
17
18module Data.Packed.Development (
19 createVector, createMatrix,
20 vec, mat,
21 app1, app2, app3, app4,
22 app5, app6, app7, app8, app9, app10,
23 MatrixOrder(..), orderOf, cmat, fmat,
24 matrixFromVector,
25 unsafeFromForeignPtr,
26 unsafeToForeignPtr,
27 check, (//),
28 at', atM', fi
29) where
30
31import Data.Packed.Internal
32
diff --git a/packages/base/src/Data/Packed/Foreign.hs b/packages/base/src/Data/Packed/Foreign.hs
deleted file mode 100644
index 1ec3694..0000000
--- a/packages/base/src/Data/Packed/Foreign.hs
+++ /dev/null
@@ -1,100 +0,0 @@
1{-# LANGUAGE MagicHash, UnboxedTuples #-}
2-- | FFI and hmatrix helpers.
3--
4-- Sample usage, to upload a perspective matrix to a shader.
5--
6-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3)
7-- @
8--
9{-# OPTIONS_HADDOCK hide #-}
10module Data.Packed.Foreign
11 ( app
12 , appVector, appVectorLen
13 , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
14 , unsafeMatrixToVector, unsafeMatrixToForeignPtr
15 ) where
16import Data.Packed.Internal
17import qualified Data.Vector.Storable as S
18import Foreign (Ptr, ForeignPtr, Storable)
19import Foreign.C.Types (CInt)
20import GHC.Base (IO(..), realWorld#)
21
22{-# INLINE unsafeInlinePerformIO #-}
23-- | If we use unsafePerformIO, it may not get inlined, so in a function that returns IO (which are all safe uses of app* in this module), there would be
24-- unecessary calls to unsafePerformIO or its internals.
25unsafeInlinePerformIO :: IO a -> a
26unsafeInlinePerformIO (IO f) = case f realWorld# of
27 (# _, x #) -> x
28
29{-# INLINE app #-}
30-- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
31-- e.g.
32--
33-- @
34-- someFunction
35-- \`appMatrixLen\` m
36-- \`appVectorLen\` v
37-- \`app\` other
38-- \`app\` arguments
39-- \`app\` go here
40-- @
41--
42-- One could also write:
43--
44-- @
45-- (someFunction
46-- \`appMatrixLen\` m
47-- \`appVectorLen\` v)
48-- other
49-- arguments
50-- (go here)
51-- @
52--
53app :: (a -> b) -> a -> b
54app f = f
55
56{-# INLINE appVector #-}
57appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
58appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
59
60{-# INLINE appVectorLen #-}
61appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
62appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
63
64{-# INLINE appMatrix #-}
65appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
66appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
67
68{-# INLINE appMatrixLen #-}
69appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
70appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
71 where
72 r = fromIntegral (rows x)
73 c = fromIntegral (cols x)
74
75{-# INLINE appMatrixRaw #-}
76appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
77appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
78
79{-# INLINE appMatrixRawLen #-}
80appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
81appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
82 where
83 r = fromIntegral (rows x)
84 c = fromIntegral (cols x)
85
86infixl 1 `app`
87infixl 1 `appVector`
88infixl 1 `appMatrix`
89infixl 1 `appMatrixRaw`
90
91{-# INLINE unsafeMatrixToVector #-}
92-- | This will disregard the order of the matrix, and simply return it as-is.
93-- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
94unsafeMatrixToVector :: Matrix a -> Vector a
95unsafeMatrixToVector = xdat
96
97{-# INLINE unsafeMatrixToForeignPtr #-}
98unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
99unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
100
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs
deleted file mode 100644
index 85f1b37..0000000
--- a/packages/base/src/Data/Packed/IO.hs
+++ /dev/null
@@ -1,167 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.IO
4-- Copyright : (c) Alberto Ruiz 2010
5-- License : BSD3
6--
7-- Maintainer : Alberto Ruiz
8-- Stability : provisional
9--
10-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix'
11--
12-----------------------------------------------------------------------------
13{-# OPTIONS_HADDOCK hide #-}
14
15module Data.Packed.IO (
16 dispf, disps, dispcf, vecdisp, latexFormat, format,
17 readMatrix, fromArray2D, loadMatrix, loadMatrix', saveMatrix
18) where
19
20import Data.Packed
21import Text.Printf(printf)
22import Data.List(intersperse)
23import Data.Complex
24import Numeric.Vectorized(vectorScan,saveMatrix)
25import Control.Applicative((<$>))
26import Data.Packed.Internal
27
28{- | Creates a string from a matrix given a separator and a function to show each entry. Using
29this function the user can easily define any desired display function:
30
31@import Text.Printf(printf)@
32
33@disp = putStr . format \" \" (printf \"%.2f\")@
34
35-}
36format :: (Element t) => String -> (t -> String) -> Matrix t -> String
37format sep f m = table sep . map (map f) . toLists $ m
38
39{- | Show a matrix with \"autoscaling\" and a given number of decimal places.
40
41>>> putStr . disps 2 $ 120 * (3><4) [1..]
423x4 E3
43 0.12 0.24 0.36 0.48
44 0.60 0.72 0.84 0.96
45 1.08 1.20 1.32 1.44
46
47-}
48disps :: Int -> Matrix Double -> String
49disps d x = sdims x ++ " " ++ formatScaled d x
50
51{- | Show a matrix with a given number of decimal places.
52
53>>> dispf 2 (1/3 + ident 3)
54"3x3\n1.33 0.33 0.33\n0.33 1.33 0.33\n0.33 0.33 1.33\n"
55
56>>> putStr . dispf 2 $ (3><4)[1,1.5..]
573x4
581.00 1.50 2.00 2.50
593.00 3.50 4.00 4.50
605.00 5.50 6.00 6.50
61
62>>> putStr . unlines . tail . lines . dispf 2 . asRow $ linspace 10 (0,1)
630.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00
64
65-}
66dispf :: Int -> Matrix Double -> String
67dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
68
69sdims x = show (rows x) ++ "x" ++ show (cols x)
70
71formatFixed d x = format " " (printf ("%."++show d++"f")) $ x
72
73isInt = all lookslikeInt . toList . flatten
74
75formatScaled dec t = "E"++show o++"\n" ++ ss
76 where ss = format " " (printf fmt. g) t
77 g x | o >= 0 = x/10^(o::Int)
78 | otherwise = x*10^(-o)
79 o | rows t == 0 || cols t == 0 = 0
80 | otherwise = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
81 fmt = '%':show (dec+3) ++ '.':show dec ++"f"
82
83{- | Show a vector using a function for showing matrices.
84
85>>> putStr . vecdisp (dispf 2) $ linspace 10 (0,1)
8610 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00
87
88-}
89vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
90vecdisp f v
91 = ((show (dim v) ++ " |> ") ++) . (++"\n")
92 . unwords . lines . tail . dropWhile (not . (`elem` " \n"))
93 . f . trans . reshape 1
94 $ v
95
96{- | Tool to display matrices with latex syntax.
97
98>>> latexFormat "bmatrix" (dispf 2 $ ident 2)
99"\\begin{bmatrix}\n1 & 0\n\\\\\n0 & 1\n\\end{bmatrix}"
100
101-}
102latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
103 -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
104 -> String
105latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"
106 where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines
107
108-- | Pretty print a complex number with at most n decimal digits.
109showComplex :: Int -> Complex Double -> String
110showComplex d (a:+b)
111 | isZero a && isZero b = "0"
112 | isZero b = sa
113 | isZero a && isOne b = s2++"i"
114 | isZero a = sb++"i"
115 | isOne b = sa++s3++"i"
116 | otherwise = sa++s1++sb++"i"
117 where
118 sa = shcr d a
119 sb = shcr d b
120 s1 = if b<0 then "" else "+"
121 s2 = if b<0 then "-" else ""
122 s3 = if b<0 then "-" else "+"
123
124shcr d a | lookslikeInt a = printf "%.0f" a
125 | otherwise = printf ("%."++show d++"f") a
126
127
128lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
129 where shx = show x
130
131isZero x = show x `elem` ["0.0","-0.0"]
132isOne x = show x `elem` ["1.0","-1.0"]
133
134-- | Pretty print a complex matrix with at most n decimal digits.
135dispcf :: Int -> Matrix (Complex Double) -> String
136dispcf d m = sdims m ++ "\n" ++ format " " (showComplex d) m
137
138--------------------------------------------------------------------
139
140-- | reads a matrix from a string containing a table of numbers.
141readMatrix :: String -> Matrix Double
142readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
143
144--------------------------------------------------------------------------------
145
146apparentCols :: FilePath -> IO Int
147apparentCols s = f . dropWhile null . map words . lines <$> readFile s
148 where
149 f [] = 0
150 f (x:_) = length x
151
152
153-- | load a matrix from an ASCII file formatted as a 2D table.
154loadMatrix :: FilePath -> IO (Matrix Double)
155loadMatrix f = do
156 v <- vectorScan f
157 c <- apparentCols f
158 if (dim v `mod` c /= 0)
159 then
160 error $ printf "loadMatrix: %d elements and %d columns in file %s"
161 (dim v) c f
162 else
163 return (reshape c v)
164
165
166loadMatrix' name = mbCatch (loadMatrix name)
167
diff --git a/packages/base/src/Data/Packed/Internal.hs b/packages/base/src/Data/Packed/Internal.hs
deleted file mode 100644
index 59a72fc..0000000
--- a/packages/base/src/Data/Packed/Internal.hs
+++ /dev/null
@@ -1,24 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Data.Packed.Internal
4-- Copyright : (c) Alberto Ruiz 2007
5-- License : BSD3
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-- Reexports all internal modules
10--
11-----------------------------------------------------------------------------
12-- #hide
13
14module Data.Packed.Internal (
15 module Data.Packed.Internal.Common,
16 module Data.Packed.Internal.Signatures,
17 module Data.Packed.Internal.Vector,
18 module Data.Packed.Internal.Matrix,
19) where
20
21import Data.Packed.Internal.Common
22import Data.Packed.Internal.Signatures
23import Data.Packed.Internal.Vector
24import Data.Packed.Internal.Matrix
diff --git a/packages/base/src/Data/Packed/Internal/Common.hs b/packages/base/src/Data/Packed/Internal/Common.hs
deleted file mode 100644
index 615bbdf..0000000
--- a/packages/base/src/Data/Packed/Internal/Common.hs
+++ /dev/null
@@ -1,160 +0,0 @@
1{-# LANGUAGE CPP #-}
2-- |
3-- Module : Data.Packed.Internal.Common
4-- Copyright : (c) Alberto Ruiz 2007
5-- License : BSD3
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9--
10-- Development utilities.
11--
12
13
14module Data.Packed.Internal.Common(
15 Adapt,
16 app1, app2, app3, app4,
17 app5, app6, app7, app8, app9, app10,
18 (//), check, mbCatch,
19 splitEvery, common, compatdim,
20 fi,
21 table,
22 finit
23) where
24
25import Control.Monad(when)
26import Foreign.C.Types
27import Foreign.Storable.Complex()
28import Data.List(transpose,intersperse)
29import Control.Exception as E
30
31-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
32splitEvery :: Int -> [a] -> [[a]]
33splitEvery _ [] = []
34splitEvery k l = take k l : splitEvery k (drop k l)
35
36-- | obtains the common value of a property of a list
37common :: (Eq a) => (b->a) -> [b] -> Maybe a
38common f = commonval . map f where
39 commonval :: (Eq a) => [a] -> Maybe a
40 commonval [] = Nothing
41 commonval [a] = Just a
42 commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
43
44-- | common value with \"adaptable\" 1
45compatdim :: [Int] -> Maybe Int
46compatdim [] = Nothing
47compatdim [a] = Just a
48compatdim (a:b:xs)
49 | a==b = compatdim (b:xs)
50 | a==1 = compatdim (b:xs)
51 | b==1 = compatdim (a:xs)
52 | otherwise = Nothing
53
54-- | Formatting tool
55table :: String -> [[String]] -> String
56table sep as = unlines . map unwords' $ transpose mtp where
57 mt = transpose as
58 longs = map (maximum . map length) mt
59 mtp = zipWith (\a b -> map (pad a) b) longs mt
60 pad n str = replicate (n - length str) ' ' ++ str
61 unwords' = concat . intersperse sep
62
63-- | postfix function application (@flip ($)@)
64(//) :: x -> (x -> y) -> y
65infixl 0 //
66(//) = flip ($)
67
68-- | specialized fromIntegral
69fi :: Int -> CInt
70fi = fromIntegral
71
72-- hmm..
73ww2 w1 o1 w2 o2 f = w1 o1 $ w2 o2 . f
74ww3 w1 o1 w2 o2 w3 o3 f = w1 o1 $ ww2 w2 o2 w3 o3 . f
75ww4 w1 o1 w2 o2 w3 o3 w4 o4 f = w1 o1 $ ww3 w2 o2 w3 o3 w4 o4 . f
76ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 f = w1 o1 $ ww4 w2 o2 w3 o3 w4 o4 w5 o5 . f
77ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 f = w1 o1 $ ww5 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 . f
78ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 f = w1 o1 $ ww6 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 . f
79ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 f = w1 o1 $ ww7 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 . f
80ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 f = w1 o1 $ ww8 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 . f
81ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 f = w1 o1 $ ww9 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 . f
82
83type Adapt f t r = t -> ((f -> r) -> IO()) -> IO()
84
85type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO()
86type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2
87type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3
88type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4
89type Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 = Adapt f t1 r1 -> t1 -> Adapt4 r1 t2 r2 t3 r3 t4 r4 t5
90type Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 = Adapt f t1 r1 -> t1 -> Adapt5 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
91type Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 = Adapt f t1 r1 -> t1 -> Adapt6 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
92type Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 = Adapt f t1 r1 -> t1 -> Adapt7 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
93type Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 = Adapt f t1 r1 -> t1 -> Adapt8 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
94type Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 = Adapt f t1 r1 -> t1 -> Adapt9 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
95
96app1 :: f -> Adapt1 f t1
97app2 :: f -> Adapt2 f t1 r1 t2
98app3 :: f -> Adapt3 f t1 r1 t2 r2 t3
99app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4
100app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5
101app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6
102app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7
103app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8
104app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9
105app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10
106
107app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s
108app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s
109app3 f w1 o1 w2 o2 w3 o3 s = ww3 w1 o1 w2 o2 w3 o3 $
110 \a1 a2 a3 -> f // a1 // a2 // a3 // check s
111app4 f w1 o1 w2 o2 w3 o3 w4 o4 s = ww4 w1 o1 w2 o2 w3 o3 w4 o4 $
112 \a1 a2 a3 a4 -> f // a1 // a2 // a3 // a4 // check s
113app5 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 s = ww5 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 $
114 \a1 a2 a3 a4 a5 -> f // a1 // a2 // a3 // a4 // a5 // check s
115app6 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 s = ww6 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 $
116 \a1 a2 a3 a4 a5 a6 -> f // a1 // a2 // a3 // a4 // a5 // a6 // check s
117app7 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 s = ww7 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 $
118 \a1 a2 a3 a4 a5 a6 a7 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // check s
119app8 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 s = ww8 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 $
120 \a1 a2 a3 a4 a5 a6 a7 a8 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // check s
121app9 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 s = ww9 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 $
122 \a1 a2 a3 a4 a5 a6 a7 a8 a9 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // check s
123app10 f w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 s = ww10 w1 o1 w2 o2 w3 o3 w4 o4 w5 o5 w6 o6 w7 o7 w8 o8 w9 o9 w10 o10 $
124 \a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 -> f // a1 // a2 // a3 // a4 // a5 // a6 // a7 // a8 // a9 // a10 // check s
125
126
127
128-- GSL error codes are <= 1024
129-- | error codes for the auxiliary functions required by the wrappers
130errorCode :: CInt -> String
131errorCode 2000 = "bad size"
132errorCode 2001 = "bad function code"
133errorCode 2002 = "memory problem"
134errorCode 2003 = "bad file"
135errorCode 2004 = "singular"
136errorCode 2005 = "didn't converge"
137errorCode 2006 = "the input matrix is not positive definite"
138errorCode 2007 = "not yet supported in this OS"
139errorCode n = "code "++show n
140
141
142-- | clear the fpu
143foreign import ccall unsafe "asm_finit" finit :: IO ()
144
145-- | check the error code
146check :: String -> IO CInt -> IO ()
147check msg f = do
148#if FINIT
149 finit
150#endif
151 err <- f
152 when (err/=0) $ error (msg++": "++errorCode err)
153 return ()
154
155-- | Error capture and conversion to Maybe
156mbCatch :: IO x -> IO (Maybe x)
157mbCatch act = E.catch (Just `fmap` act) f
158 where f :: SomeException -> IO (Maybe x)
159 f _ = return Nothing
160
diff --git a/packages/base/src/Data/Packed/Internal/Matrix.hs b/packages/base/src/Data/Packed/Internal/Matrix.hs
deleted file mode 100644
index 150b978..0000000
--- a/packages/base/src/Data/Packed/Internal/Matrix.hs
+++ /dev/null
@@ -1,423 +0,0 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE BangPatterns #-}
5
6-- |
7-- Module : Data.Packed.Internal.Matrix
8-- Copyright : (c) Alberto Ruiz 2007
9-- License : BSD3
10-- Maintainer : Alberto Ruiz
11-- Stability : provisional
12--
13-- Internal matrix representation
14--
15
16module Data.Packed.Internal.Matrix(
17 Matrix(..), rows, cols, cdat, fdat,
18 MatrixOrder(..), orderOf,
19 createMatrix, mat,
20 cmat, fmat,
21 toLists, flatten, reshape,
22 Element(..),
23 trans,
24 fromRows, toRows, fromColumns, toColumns,
25 matrixFromVector,
26 subMatrix,
27 liftMatrix, liftMatrix2,
28 (@@>), atM',
29 singleton,
30 emptyM,
31 size, shSize, conformVs, conformMs, conformVTo, conformMTo
32) where
33
34import Data.Packed.Internal.Common
35import Data.Packed.Internal.Signatures
36import Data.Packed.Internal.Vector
37
38import Foreign.Marshal.Alloc(alloca, free)
39import Foreign.Marshal.Array(newArray)
40import Foreign.Ptr(Ptr, castPtr)
41import Foreign.Storable(Storable, peekElemOff, pokeElemOff, poke, sizeOf)
42import Data.Complex(Complex)
43import Foreign.C.Types
44import System.IO.Unsafe(unsafePerformIO)
45import Control.DeepSeq
46
47-----------------------------------------------------------------
48
49{- Design considerations for the Matrix Type
50 -----------------------------------------
51
52- we must easily handle both row major and column major order,
53 for bindings to LAPACK and GSL/C
54
55- we'd like to simplify redundant matrix transposes:
56 - Some of them arise from the order requirements of some functions
57 - some functions (matrix product) admit transposed arguments
58
59- maybe we don't really need this kind of simplification:
60 - more complex code
61 - some computational overhead
62 - only appreciable gain in code with a lot of redundant transpositions
63 and cheap matrix computations
64
65- we could carry both the matrix and its (lazily computed) transpose.
66 This may save some transpositions, but it is necessary to keep track of the
67 data which is actually computed to be used by functions like the matrix product
68 which admit both orders.
69
70- but if we need the transposed data and it is not in the structure, we must make
71 sure that we touch the same foreignptr that is used in the computation.
72
73- a reasonable solution is using two constructors for a matrix. Transposition just
74 "flips" the constructor. Actual data transposition is not done if followed by a
75 matrix product or another transpose.
76
77-}
78
79data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
80
81transOrder RowMajor = ColumnMajor
82transOrder ColumnMajor = RowMajor
83{- | Matrix representation suitable for BLAS\/LAPACK computations.
84
85The elements are stored in a continuous memory array.
86
87-}
88
89data Matrix t = Matrix { irows :: {-# UNPACK #-} !Int
90 , icols :: {-# UNPACK #-} !Int
91 , xdat :: {-# UNPACK #-} !(Vector t)
92 , order :: !MatrixOrder }
93-- RowMajor: preferred by C, fdat may require a transposition
94-- ColumnMajor: preferred by LAPACK, cdat may require a transposition
95
96cdat = xdat
97fdat = xdat
98
99rows :: Matrix t -> Int
100rows = irows
101
102cols :: Matrix t -> Int
103cols = icols
104
105orderOf :: Matrix t -> MatrixOrder
106orderOf = order
107
108
109-- | Matrix transpose.
110trans :: Matrix t -> Matrix t
111trans Matrix {irows = r, icols = c, xdat = d, order = o } = Matrix { irows = c, icols = r, xdat = d, order = transOrder o}
112
113cmat :: (Element t) => Matrix t -> Matrix t
114cmat m@Matrix{order = RowMajor} = m
115cmat Matrix {irows = r, icols = c, xdat = d, order = ColumnMajor } = Matrix { irows = r, icols = c, xdat = transdata r d c, order = RowMajor}
116
117fmat :: (Element t) => Matrix t -> Matrix t
118fmat m@Matrix{order = ColumnMajor} = m
119fmat Matrix {irows = r, icols = c, xdat = d, order = RowMajor } = Matrix { irows = r, icols = c, xdat = transdata c d r, order = ColumnMajor}
120
121-- C-Haskell matrix adapter
122-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
123
124mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
125mat a f =
126 unsafeWith (xdat a) $ \p -> do
127 let m g = do
128 g (fi (rows a)) (fi (cols a)) p
129 f m
130
131{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
132
133>>> flatten (ident 3)
134fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
135
136-}
137flatten :: Element t => Matrix t -> Vector t
138flatten = xdat . cmat
139
140{-
141type Mt t s = Int -> Int -> Ptr t -> s
142
143infixr 6 ::>
144type t ::> s = Mt t s
145-}
146
147-- | the inverse of 'Data.Packed.Matrix.fromLists'
148toLists :: (Element t) => Matrix t -> [[t]]
149toLists m = splitEvery (cols m) . toList . flatten $ m
150
151-- | Create a matrix from a list of vectors.
152-- All vectors must have the same dimension,
153-- or dimension 1, which is are automatically expanded.
154fromRows :: Element t => [Vector t] -> Matrix t
155fromRows [] = emptyM 0 0
156fromRows vs = case compatdim (map dim vs) of
157 Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs)
158 Just 0 -> emptyM r 0
159 Just c -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs
160 where
161 r = length vs
162 adapt c v
163 | c == 0 = fromList[]
164 | dim v == c = v
165 | otherwise = constantD (v@>0) c
166
167-- | extracts the rows of a matrix as a list of vectors
168toRows :: Element t => Matrix t -> [Vector t]
169toRows m
170 | c == 0 = replicate r (fromList[])
171 | otherwise = toRows' 0
172 where
173 v = flatten m
174 r = rows m
175 c = cols m
176 toRows' k | k == r*c = []
177 | otherwise = subVector k c v : toRows' (k+c)
178
179-- | Creates a matrix from a list of vectors, as columns
180fromColumns :: Element t => [Vector t] -> Matrix t
181fromColumns m = trans . fromRows $ m
182
183-- | Creates a list of vectors from the columns of a matrix
184toColumns :: Element t => Matrix t -> [Vector t]
185toColumns m = toRows . trans $ m
186
187-- | Reads a matrix position.
188(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
189infixl 9 @@>
190m@Matrix {irows = r, icols = c} @@> (i,j)
191 | safe = if i<0 || i>=r || j<0 || j>=c
192 then error "matrix indexing out of range"
193 else atM' m i j
194 | otherwise = atM' m i j
195{-# INLINE (@@>) #-}
196
197-- Unsafe matrix access without range checking
198atM' Matrix {icols = c, xdat = v, order = RowMajor} i j = v `at'` (i*c+j)
199atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i)
200{-# INLINE atM' #-}
201
202------------------------------------------------------------------
203
204matrixFromVector o r c v
205 | r * c == dim v = m
206 | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
207 where
208 m = Matrix { irows = r, icols = c, xdat = v, order = o }
209
210-- allocates memory for a new matrix
211createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a)
212createMatrix ord r c = do
213 p <- createVector (r*c)
214 return (matrixFromVector ord r c p)
215
216{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@
217where r is the desired number of rows.)
218
219>>> reshape 4 (fromList [1..12])
220(3><4)
221 [ 1.0, 2.0, 3.0, 4.0
222 , 5.0, 6.0, 7.0, 8.0
223 , 9.0, 10.0, 11.0, 12.0 ]
224
225-}
226reshape :: Storable t => Int -> Vector t -> Matrix t
227reshape 0 v = matrixFromVector RowMajor 0 0 v
228reshape c v = matrixFromVector RowMajor (dim v `div` c) c v
229
230singleton x = reshape 1 (fromList [x])
231
232-- | application of a vector function on the flattened matrix elements
233liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
234liftMatrix f Matrix { irows = r, icols = c, xdat = d, order = o } = matrixFromVector o r c (f d)
235
236-- | application of a vector function on the flattened matrices elements
237liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
238liftMatrix2 f m1 m2
239 | not (compat m1 m2) = error "nonconformant matrices in liftMatrix2"
240 | otherwise = case orderOf m1 of
241 RowMajor -> matrixFromVector RowMajor (rows m1) (cols m1) (f (xdat m1) (flatten m2))
242 ColumnMajor -> matrixFromVector ColumnMajor (rows m1) (cols m1) (f (xdat m1) ((xdat.fmat) m2))
243
244
245compat :: Matrix a -> Matrix b -> Bool
246compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
247
248------------------------------------------------------------------
249
250{- | Supported matrix elements.
251
252 This class provides optimized internal
253 operations for selected element types.
254 It provides unoptimised defaults for any 'Storable' type,
255 so you can create instances simply as:
256
257 >instance Element Foo
258-}
259class (Storable a) => Element a where
260 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position
261 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
262 -> Matrix a -> Matrix a
263 subMatrixD = subMatrix'
264 transdata :: Int -> Vector a -> Int -> Vector a
265 transdata = transdataP -- transdata'
266 constantD :: a -> Int -> Vector a
267 constantD = constantP -- constant'
268
269
270instance Element Float where
271 transdata = transdataAux ctransF
272 constantD = constantAux cconstantF
273
274instance Element Double where
275 transdata = transdataAux ctransR
276 constantD = constantAux cconstantR
277
278instance Element (Complex Float) where
279 transdata = transdataAux ctransQ
280 constantD = constantAux cconstantQ
281
282instance Element (Complex Double) where
283 transdata = transdataAux ctransC
284 constantD = constantAux cconstantC
285
286-------------------------------------------------------------------
287
288transdataAux fun c1 d c2 =
289 if noneed
290 then d
291 else unsafePerformIO $ do
292 v <- createVector (dim d)
293 unsafeWith d $ \pd ->
294 unsafeWith v $ \pv ->
295 fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
296 return v
297 where r1 = dim d `div` c1
298 r2 = dim d `div` c2
299 noneed = dim d == 0 || r1 == 1 || c1 == 1
300
301transdataP :: Storable a => Int -> Vector a -> Int -> Vector a
302transdataP c1 d c2 =
303 if noneed
304 then d
305 else unsafePerformIO $ do
306 v <- createVector (dim d)
307 unsafeWith d $ \pd ->
308 unsafeWith v $ \pv ->
309 ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP"
310 return v
311 where r1 = dim d `div` c1
312 r2 = dim d `div` c2
313 sz = sizeOf (d @> 0)
314 noneed = dim d == 0 || r1 == 1 || c1 == 1
315
316foreign import ccall unsafe "transF" ctransF :: TFMFM
317foreign import ccall unsafe "transR" ctransR :: TMM
318foreign import ccall unsafe "transQ" ctransQ :: TQMQM
319foreign import ccall unsafe "transC" ctransC :: TCMCM
320foreign import ccall unsafe "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt
321
322----------------------------------------------------------------------
323
324constantAux fun x n = unsafePerformIO $ do
325 v <- createVector n
326 px <- newArray [x]
327 app1 (fun px) vec v "constantAux"
328 free px
329 return v
330
331foreign import ccall unsafe "constantF" cconstantF :: Ptr Float -> TF
332
333foreign import ccall unsafe "constantR" cconstantR :: Ptr Double -> TV
334
335foreign import ccall unsafe "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV
336
337foreign import ccall unsafe "constantC" cconstantC :: Ptr (Complex Double) -> TCV
338
339constantP :: Storable a => a -> Int -> Vector a
340constantP a n = unsafePerformIO $ do
341 let sz = sizeOf a
342 v <- createVector n
343 unsafeWith v $ \p -> do
344 alloca $ \k -> do
345 poke k a
346 cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP"
347 return v
348foreign import ccall unsafe "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt
349
350----------------------------------------------------------------------
351
352-- | Extracts a submatrix from a matrix.
353subMatrix :: Element a
354 => (Int,Int) -- ^ (r0,c0) starting position
355 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
356 -> Matrix a -- ^ input matrix
357 -> Matrix a -- ^ result
358subMatrix (r0,c0) (rt,ct) m
359 | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) &&
360 0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
361 | otherwise = error $ "wrong subMatrix "++
362 show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
363
364subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do
365 w <- createVector (rt*ct)
366 unsafeWith v $ \p ->
367 unsafeWith w $ \q -> do
368 let go (-1) _ = return ()
369 go !i (-1) = go (i-1) (ct-1)
370 go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0)
371 pokeElemOff q (i*ct+j) x
372 go i (j-1)
373 go (rt-1) (ct-1)
374 return w
375
376subMatrix' (r0,c0) (rt,ct) (Matrix { icols = c, xdat = v, order = RowMajor}) = Matrix rt ct (subMatrix'' (r0,c0) (rt,ct) c v) RowMajor
377subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m)
378
379--------------------------------------------------------------------------
380
381maxZ xs = if minimum xs == 0 then 0 else maximum xs
382
383conformMs ms = map (conformMTo (r,c)) ms
384 where
385 r = maxZ (map rows ms)
386 c = maxZ (map cols ms)
387
388
389conformVs vs = map (conformVTo n) vs
390 where
391 n = maxZ (map dim vs)
392
393conformMTo (r,c) m
394 | size m == (r,c) = m
395 | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
396 | size m == (r,1) = repCols c m
397 | size m == (1,c) = repRows r m
398 | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
399
400conformVTo n v
401 | dim v == n = v
402 | dim v == 1 = constantD (v@>0) n
403 | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n
404
405repRows n x = fromRows (replicate n (flatten x))
406repCols n x = fromColumns (replicate n (flatten x))
407
408size m = (rows m, cols m)
409
410shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
411
412emptyM r c = matrixFromVector RowMajor r c (fromList[])
413
414----------------------------------------------------------------------
415
416instance (Storable t, NFData t) => NFData (Matrix t)
417 where
418 rnf m | d > 0 = rnf (v @> 0)
419 | otherwise = ()
420 where
421 d = dim v
422 v = xdat m
423
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
deleted file mode 100644
index 257ad73..0000000
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ /dev/null
@@ -1,720 +0,0 @@
1{-# LANGUAGE CPP #-}
2{-# LANGUAGE TypeFamilies #-}
3{-# LANGUAGE FlexibleContexts #-}
4{-# LANGUAGE FlexibleInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-}
6{-# LANGUAGE FunctionalDependencies #-}
7{-# LANGUAGE UndecidableInstances #-}
8
9-----------------------------------------------------------------------------
10-- |
11-- Module : Data.Packed.Internal.Numeric
12-- Copyright : (c) Alberto Ruiz 2010-14
13-- License : BSD3
14-- Maintainer : Alberto Ruiz
15-- Stability : provisional
16--
17-----------------------------------------------------------------------------
18
19module Data.Packed.Internal.Numeric (
20 -- * Basic functions
21 ident, diag, ctrans,
22 -- * Generic operations
23 Container(..),
24 scalar, conj, scale, arctan2, cmap,
25 atIndex, minIndex, maxIndex, minElement, maxElement,
26 sumElements, prodElements,
27 step, cond, find, assoc, accum,
28 Transposable(..), Linear(..), Testable(..),
29 -- * Matrix product and related functions
30 Product(..), udot,
31 mXm,mXv,vXm,
32 outer, kronecker,
33 -- * sorting
34 sortVector,
35 -- * Element conversion
36 Convert(..),
37 Complexable(),
38 RealElement(),
39 roundVector,
40 RealOf, ComplexOf, SingleOf, DoubleOf,
41 IndexOf,
42 module Data.Complex
43) where
44
45import Data.Packed
46import Data.Packed.ST as ST
47import Numeric.Conversion
48import Data.Packed.Development
49import Numeric.Vectorized
50import Data.Complex
51import Control.Applicative((<*>))
52
53import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
54import Data.Packed.Internal
55
56-------------------------------------------------------------------
57
58type family IndexOf (c :: * -> *)
59
60type instance IndexOf Vector = Int
61type instance IndexOf Matrix = (Int,Int)
62
63type family ArgOf (c :: * -> *) a
64
65type instance ArgOf Vector a = a -> a
66type instance ArgOf Matrix a = a -> a -> a
67
68-------------------------------------------------------------------
69
70-- | Basic element-by-element functions for numeric containers
71class (Complexable c, Fractional e, Element e) => Container c e
72 where
73 size' :: c e -> IndexOf c
74 scalar' :: e -> c e
75 conj' :: c e -> c e
76 scale' :: e -> c e -> c e
77 -- | scale the element by element reciprocal of the object:
78 --
79 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
80 scaleRecip :: e -> c e -> c e
81 addConstant :: e -> c e -> c e
82 add :: c e -> c e -> c e
83 sub :: c e -> c e -> c e
84 -- | element by element multiplication
85 mul :: c e -> c e -> c e
86 -- | element by element division
87 divide :: c e -> c e -> c e
88 equal :: c e -> c e -> Bool
89 --
90 -- element by element inverse tangent
91 arctan2' :: c e -> c e -> c e
92 cmap' :: (Element b) => (e -> b) -> c e -> c b
93 konst' :: e -> IndexOf c -> c e
94 build' :: IndexOf c -> (ArgOf c e) -> c e
95 atIndex' :: c e -> IndexOf c -> e
96 minIndex' :: c e -> IndexOf c
97 maxIndex' :: c e -> IndexOf c
98 minElement' :: c e -> e
99 maxElement' :: c e -> e
100 sumElements' :: c e -> e
101 prodElements' :: c e -> e
102 step' :: RealElement e => c e -> c e
103 cond' :: RealElement e
104 => c e -- ^ a
105 -> c e -- ^ b
106 -> c e -- ^ l
107 -> c e -- ^ e
108 -> c e -- ^ g
109 -> c e -- ^ result
110 find' :: (e -> Bool) -> c e -> [IndexOf c]
111 assoc' :: IndexOf c -- ^ size
112 -> e -- ^ default value
113 -> [(IndexOf c, e)] -- ^ association list
114 -> c e -- ^ result
115 accum' :: c e -- ^ initial structure
116 -> (e -> e -> e) -- ^ update function
117 -> [(IndexOf c, e)] -- ^ association list
118 -> c e -- ^ result
119
120--------------------------------------------------------------------------
121
122instance Container Vector Float
123 where
124 size' = dim
125 scale' = vectorMapValF Scale
126 scaleRecip = vectorMapValF Recip
127 addConstant = vectorMapValF AddConstant
128 add = vectorZipF Add
129 sub = vectorZipF Sub
130 mul = vectorZipF Mul
131 divide = vectorZipF Div
132 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
133 arctan2' = vectorZipF ATan2
134 scalar' x = fromList [x]
135 konst' = constantD
136 build' = buildV
137 conj' = id
138 cmap' = mapVector
139 atIndex' = (@>)
140 minIndex' = emptyErrorV "minIndex" (round . toScalarF MinIdx)
141 maxIndex' = emptyErrorV "maxIndex" (round . toScalarF MaxIdx)
142 minElement' = emptyErrorV "minElement" (toScalarF Min)
143 maxElement' = emptyErrorV "maxElement" (toScalarF Max)
144 sumElements' = sumF
145 prodElements' = prodF
146 step' = stepF
147 find' = findV
148 assoc' = assocV
149 accum' = accumV
150 cond' = condV condF
151
152instance Container Vector Double
153 where
154 size' = dim
155 scale' = vectorMapValR Scale
156 scaleRecip = vectorMapValR Recip
157 addConstant = vectorMapValR AddConstant
158 add = vectorZipR Add
159 sub = vectorZipR Sub
160 mul = vectorZipR Mul
161 divide = vectorZipR Div
162 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
163 arctan2' = vectorZipR ATan2
164 scalar' x = fromList [x]
165 konst' = constantD
166 build' = buildV
167 conj' = id
168 cmap' = mapVector
169 atIndex' = (@>)
170 minIndex' = emptyErrorV "minIndex" (round . toScalarR MinIdx)
171 maxIndex' = emptyErrorV "maxIndex" (round . toScalarR MaxIdx)
172 minElement' = emptyErrorV "minElement" (toScalarR Min)
173 maxElement' = emptyErrorV "maxElement" (toScalarR Max)
174 sumElements' = sumR
175 prodElements' = prodR
176 step' = stepD
177 find' = findV
178 assoc' = assocV
179 accum' = accumV
180 cond' = condV condD
181
182instance Container Vector (Complex Double)
183 where
184 size' = dim
185 scale' = vectorMapValC Scale
186 scaleRecip = vectorMapValC Recip
187 addConstant = vectorMapValC AddConstant
188 add = vectorZipC Add
189 sub = vectorZipC Sub
190 mul = vectorZipC Mul
191 divide = vectorZipC Div
192 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
193 arctan2' = vectorZipC ATan2
194 scalar' x = fromList [x]
195 konst' = constantD
196 build' = buildV
197 conj' = conjugateC
198 cmap' = mapVector
199 atIndex' = (@>)
200 minIndex' = emptyErrorV "minIndex" (minIndex' . fst . fromComplex . (mul <*> conj'))
201 maxIndex' = emptyErrorV "maxIndex" (maxIndex' . fst . fromComplex . (mul <*> conj'))
202 minElement' = emptyErrorV "minElement" (atIndex' <*> minIndex')
203 maxElement' = emptyErrorV "maxElement" (atIndex' <*> maxIndex')
204 sumElements' = sumC
205 prodElements' = prodC
206 step' = undefined -- cannot match
207 find' = findV
208 assoc' = assocV
209 accum' = accumV
210 cond' = undefined -- cannot match
211
212instance Container Vector (Complex Float)
213 where
214 size' = dim
215 scale' = vectorMapValQ Scale
216 scaleRecip = vectorMapValQ Recip
217 addConstant = vectorMapValQ AddConstant
218 add = vectorZipQ Add
219 sub = vectorZipQ Sub
220 mul = vectorZipQ Mul
221 divide = vectorZipQ Div
222 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
223 arctan2' = vectorZipQ ATan2
224 scalar' x = fromList [x]
225 konst' = constantD
226 build' = buildV
227 conj' = conjugateQ
228 cmap' = mapVector
229 atIndex' = (@>)
230 minIndex' = emptyErrorV "minIndex" (minIndex' . fst . fromComplex . (mul <*> conj'))
231 maxIndex' = emptyErrorV "maxIndex" (maxIndex' . fst . fromComplex . (mul <*> conj'))
232 minElement' = emptyErrorV "minElement" (atIndex' <*> minIndex')
233 maxElement' = emptyErrorV "maxElement" (atIndex' <*> maxIndex')
234 sumElements' = sumQ
235 prodElements' = prodQ
236 step' = undefined -- cannot match
237 find' = findV
238 assoc' = assocV
239 accum' = accumV
240 cond' = undefined -- cannot match
241
242---------------------------------------------------------------
243
244instance (Fractional a, Element a, Container Vector a) => Container Matrix a
245 where
246 size' = size
247 scale' x = liftMatrix (scale' x)
248 scaleRecip x = liftMatrix (scaleRecip x)
249 addConstant x = liftMatrix (addConstant x)
250 add = liftMatrix2 add
251 sub = liftMatrix2 sub
252 mul = liftMatrix2 mul
253 divide = liftMatrix2 divide
254 equal a b = cols a == cols b && flatten a `equal` flatten b
255 arctan2' = liftMatrix2 arctan2'
256 scalar' x = (1><1) [x]
257 konst' v (r,c) = matrixFromVector RowMajor r c (konst' v (r*c))
258 build' = buildM
259 conj' = liftMatrix conj'
260 cmap' f = liftMatrix (mapVector f)
261 atIndex' = (@@>)
262 minIndex' = emptyErrorM "minIndex of Matrix" $
263 \m -> divMod (minIndex' $ flatten m) (cols m)
264 maxIndex' = emptyErrorM "maxIndex of Matrix" $
265 \m -> divMod (maxIndex' $ flatten m) (cols m)
266 minElement' = emptyErrorM "minElement of Matrix" (atIndex' <*> minIndex')
267 maxElement' = emptyErrorM "maxElement of Matrix" (atIndex' <*> maxIndex')
268 sumElements' = sumElements . flatten
269 prodElements' = prodElements . flatten
270 step' = liftMatrix step
271 find' = findM
272 assoc' = assocM
273 accum' = accumM
274 cond' = condM
275
276
277emptyErrorV msg f v =
278 if dim v > 0
279 then f v
280 else error $ msg ++ " of Vector with dim = 0"
281
282emptyErrorM msg f m =
283 if rows m > 0 && cols m > 0
284 then f m
285 else error $ msg++" "++shSize m
286
287--------------------------------------------------------------------------------
288
289-- | create a structure with a single element
290--
291-- >>> let v = fromList [1..3::Double]
292-- >>> v / scalar (norm2 v)
293-- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732]
294--
295scalar :: Container c e => e -> c e
296scalar = scalar'
297
298-- | complex conjugate
299conj :: Container c e => c e -> c e
300conj = conj'
301
302-- | multiplication by scalar
303scale :: Container c e => e -> c e -> c e
304scale = scale'
305
306arctan2 :: Container c e => c e -> c e -> c e
307arctan2 = arctan2'
308
309-- | like 'fmap' (cannot implement instance Functor because of Element class constraint)
310cmap :: (Element b, Container c e) => (e -> b) -> c e -> c b
311cmap = cmap'
312
313-- | indexing function
314atIndex :: Container c e => c e -> IndexOf c -> e
315atIndex = atIndex'
316
317-- | index of minimum element
318minIndex :: Container c e => c e -> IndexOf c
319minIndex = minIndex'
320
321-- | index of maximum element
322maxIndex :: Container c e => c e -> IndexOf c
323maxIndex = maxIndex'
324
325-- | value of minimum element
326minElement :: Container c e => c e -> e
327minElement = minElement'
328
329-- | value of maximum element
330maxElement :: Container c e => c e -> e
331maxElement = maxElement'
332
333-- | the sum of elements
334sumElements :: Container c e => c e -> e
335sumElements = sumElements'
336
337-- | the product of elements
338prodElements :: Container c e => c e -> e
339prodElements = prodElements'
340
341
342-- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@
343--
344-- >>> step $ linspace 5 (-1,1::Double)
345-- 5 |> [0.0,0.0,0.0,1.0,1.0]
346--
347step
348 :: (RealElement e, Container c e)
349 => c e
350 -> c e
351step = step'
352
353
354-- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@.
355--
356-- Arguments with any dimension = 1 are automatically expanded:
357--
358-- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double
359-- (3><4)
360-- [ 100.0, 2.0, 3.0, 4.0
361-- , 0.0, 100.0, 7.0, 8.0
362-- , 0.0, 0.0, 100.0, 12.0 ]
363--
364cond
365 :: (RealElement e, Container c e)
366 => c e -- ^ a
367 -> c e -- ^ b
368 -> c e -- ^ l
369 -> c e -- ^ e
370 -> c e -- ^ g
371 -> c e -- ^ result
372cond = cond'
373
374
375-- | Find index of elements which satisfy a predicate
376--
377-- >>> find (>0) (ident 3 :: Matrix Double)
378-- [(0,0),(1,1),(2,2)]
379--
380find
381 :: Container c e
382 => (e -> Bool)
383 -> c e
384 -> [IndexOf c]
385find = find'
386
387
388-- | Create a structure from an association list
389--
390-- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double
391-- fromList [0.0,4.0,0.0,7.0,0.0]
392--
393-- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double)
394-- (2><3)
395-- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0
396-- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ]
397--
398assoc
399 :: Container c e
400 => IndexOf c -- ^ size
401 -> e -- ^ default value
402 -> [(IndexOf c, e)] -- ^ association list
403 -> c e -- ^ result
404assoc = assoc'
405
406
407-- | Modify a structure using an update function
408--
409-- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double
410-- (5><5)
411-- [ 1.0, 0.0, 0.0, 3.0, 0.0
412-- , 0.0, 6.0, 0.0, 0.0, 0.0
413-- , 0.0, 0.0, 1.0, 0.0, 0.0
414-- , 0.0, 0.0, 0.0, 1.0, 0.0
415-- , 0.0, 0.0, 0.0, 0.0, 1.0 ]
416--
417-- computation of histogram:
418--
419-- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double
420-- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0]
421--
422accum
423 :: Container c e
424 => c e -- ^ initial structure
425 -> (e -> e -> e) -- ^ update function
426 -> [(IndexOf c, e)] -- ^ association list
427 -> c e -- ^ result
428accum = accum'
429
430
431--------------------------------------------------------------------------------
432
433-- | Matrix product and related functions
434class (Num e, Element e) => Product e where
435 -- | matrix product
436 multiply :: Matrix e -> Matrix e -> Matrix e
437 -- | sum of absolute value of elements (differs in complex case from @norm1@)
438 absSum :: Vector e -> RealOf e
439 -- | sum of absolute value of elements
440 norm1 :: Vector e -> RealOf e
441 -- | euclidean norm
442 norm2 :: Vector e -> RealOf e
443 -- | element of maximum magnitude
444 normInf :: Vector e -> RealOf e
445
446instance Product Float where
447 norm2 = emptyVal (toScalarF Norm2)
448 absSum = emptyVal (toScalarF AbsSum)
449 norm1 = emptyVal (toScalarF AbsSum)
450 normInf = emptyVal (maxElement . vectorMapF Abs)
451 multiply = emptyMul multiplyF
452
453instance Product Double where
454 norm2 = emptyVal (toScalarR Norm2)
455 absSum = emptyVal (toScalarR AbsSum)
456 norm1 = emptyVal (toScalarR AbsSum)
457 normInf = emptyVal (maxElement . vectorMapR Abs)
458 multiply = emptyMul multiplyR
459
460instance Product (Complex Float) where
461 norm2 = emptyVal (toScalarQ Norm2)
462 absSum = emptyVal (toScalarQ AbsSum)
463 norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapQ Abs)
464 normInf = emptyVal (maxElement . fst . fromComplex . vectorMapQ Abs)
465 multiply = emptyMul multiplyQ
466
467instance Product (Complex Double) where
468 norm2 = emptyVal (toScalarC Norm2)
469 absSum = emptyVal (toScalarC AbsSum)
470 norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapC Abs)
471 normInf = emptyVal (maxElement . fst . fromComplex . vectorMapC Abs)
472 multiply = emptyMul multiplyC
473
474emptyMul m a b
475 | x1 == 0 && x2 == 0 || r == 0 || c == 0 = konst' 0 (r,c)
476 | otherwise = m a b
477 where
478 r = rows a
479 x1 = cols a
480 x2 = rows b
481 c = cols b
482
483emptyVal f v =
484 if dim v > 0
485 then f v
486 else 0
487
488-- FIXME remove unused C wrappers
489-- | unconjugated dot product
490udot :: Product e => Vector e -> Vector e -> e
491udot u v
492 | dim u == dim v = val (asRow u `multiply` asColumn v)
493 | otherwise = error $ "different dimensions "++show (dim u)++" and "++show (dim v)++" in dot product"
494 where
495 val m | dim u > 0 = m@@>(0,0)
496 | otherwise = 0
497
498----------------------------------------------------------
499
500-- synonym for matrix product
501mXm :: Product t => Matrix t -> Matrix t -> Matrix t
502mXm = multiply
503
504-- matrix - vector product
505mXv :: Product t => Matrix t -> Vector t -> Vector t
506mXv m v = flatten $ m `mXm` (asColumn v)
507
508-- vector - matrix product
509vXm :: Product t => Vector t -> Matrix t -> Vector t
510vXm v m = flatten $ (asRow v) `mXm` m
511
512{- | Outer product of two vectors.
513
514>>> fromList [1,2,3] `outer` fromList [5,2,3]
515(3><3)
516 [ 5.0, 2.0, 3.0
517 , 10.0, 4.0, 6.0
518 , 15.0, 6.0, 9.0 ]
519
520-}
521outer :: (Product t) => Vector t -> Vector t -> Matrix t
522outer u v = asColumn u `multiply` asRow v
523
524{- | Kronecker product of two matrices.
525
526@m1=(2><3)
527 [ 1.0, 2.0, 0.0
528 , 0.0, -1.0, 3.0 ]
529m2=(4><3)
530 [ 1.0, 2.0, 3.0
531 , 4.0, 5.0, 6.0
532 , 7.0, 8.0, 9.0
533 , 10.0, 11.0, 12.0 ]@
534
535>>> kronecker m1 m2
536(8><9)
537 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
538 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
539 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
540 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
541 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
542 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
543 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
544 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]
545
546-}
547kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
548kronecker a b = fromBlocks
549 . splitEvery (cols a)
550 . map (reshape (cols b))
551 . toRows
552 $ flatten a `outer` flatten b
553
554-------------------------------------------------------------------
555
556
557class Convert t where
558 real :: Container c t => c (RealOf t) -> c t
559 complex :: Container c t => c t -> c (ComplexOf t)
560 single :: Container c t => c t -> c (SingleOf t)
561 double :: Container c t => c t -> c (DoubleOf t)
562 toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)
563 fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)
564
565
566instance Convert Double where
567 real = id
568 complex = comp'
569 single = single'
570 double = id
571 toComplex = toComplex'
572 fromComplex = fromComplex'
573
574instance Convert Float where
575 real = id
576 complex = comp'
577 single = id
578 double = double'
579 toComplex = toComplex'
580 fromComplex = fromComplex'
581
582instance Convert (Complex Double) where
583 real = comp'
584 complex = id
585 single = single'
586 double = id
587 toComplex = toComplex'
588 fromComplex = fromComplex'
589
590instance Convert (Complex Float) where
591 real = comp'
592 complex = id
593 single = id
594 double = double'
595 toComplex = toComplex'
596 fromComplex = fromComplex'
597
598-------------------------------------------------------------------
599
600type family RealOf x
601
602type instance RealOf Double = Double
603type instance RealOf (Complex Double) = Double
604
605type instance RealOf Float = Float
606type instance RealOf (Complex Float) = Float
607
608type family ComplexOf x
609
610type instance ComplexOf Double = Complex Double
611type instance ComplexOf (Complex Double) = Complex Double
612
613type instance ComplexOf Float = Complex Float
614type instance ComplexOf (Complex Float) = Complex Float
615
616type family SingleOf x
617
618type instance SingleOf Double = Float
619type instance SingleOf Float = Float
620
621type instance SingleOf (Complex a) = Complex (SingleOf a)
622
623type family DoubleOf x
624
625type instance DoubleOf Double = Double
626type instance DoubleOf Float = Double
627
628type instance DoubleOf (Complex a) = Complex (DoubleOf a)
629
630type family ElementOf c
631
632type instance ElementOf (Vector a) = a
633type instance ElementOf (Matrix a) = a
634
635------------------------------------------------------------
636
637buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]
638 where rs = map fromIntegral [0 .. (rc-1)]
639 cs = map fromIntegral [0 .. (cc-1)]
640
641buildV n f = fromList [f k | k <- ks]
642 where ks = map fromIntegral [0 .. (n-1)]
643
644--------------------------------------------------------
645-- | conjugate transpose
646ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e
647ctrans = liftMatrix conj' . trans
648
649-- | Creates a square matrix with a given diagonal.
650diag :: (Num a, Element a) => Vector a -> Matrix a
651diag v = diagRect 0 v n n where n = dim v
652
653-- | creates the identity matrix of given dimension
654ident :: (Num a, Element a) => Int -> Matrix a
655ident n = diag (constantD 1 n)
656
657--------------------------------------------------------
658
659findV p x = foldVectorWithIndex g [] x where
660 g k z l = if p z then k:l else l
661
662findM p x = map ((`divMod` cols x)) $ findV p (flatten x)
663
664assocV n z xs = ST.runSTVector $ do
665 v <- ST.newVector z n
666 mapM_ (\(k,x) -> ST.writeVector v k x) xs
667 return v
668
669assocM (r,c) z xs = ST.runSTMatrix $ do
670 m <- ST.newMatrix z r c
671 mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs
672 return m
673
674accumV v0 f xs = ST.runSTVector $ do
675 v <- ST.thawVector v0
676 mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs
677 return v
678
679accumM m0 f xs = ST.runSTMatrix $ do
680 m <- ST.thawMatrix m0
681 mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs
682 return m
683
684----------------------------------------------------------------------
685
686condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t'
687 where
688 args@(a'':_) = conformMs [a,b,l,e,t]
689 [a', b', l', e', t'] = map flatten args
690
691condV f a b l e t = f a' b' l' e' t'
692 where
693 [a', b', l', e', t'] = conformVs [a,b,l,e,t]
694
695--------------------------------------------------------------------------------
696
697class Transposable m mt | m -> mt, mt -> m
698 where
699 -- | (conjugate) transpose
700 tr :: m -> mt
701
702instance (Container Vector t) => Transposable (Matrix t) (Matrix t)
703 where
704 tr = ctrans
705
706class Linear t v
707 where
708 scalarL :: t -> v
709 addL :: v -> v -> v
710 scaleL :: t -> v -> v
711
712
713class Testable t
714 where
715 checkT :: t -> (Bool, IO())
716 ioCheckT :: t -> IO (Bool, IO())
717 ioCheckT = return . checkT
718
719--------------------------------------------------------------------------------
720
diff --git a/packages/base/src/Data/Packed/Internal/Signatures.hs b/packages/base/src/Data/Packed/Internal/Signatures.hs
deleted file mode 100644
index acc3070..0000000
--- a/packages/base/src/Data/Packed/Internal/Signatures.hs
+++ /dev/null
@@ -1,70 +0,0 @@
1-- |
2-- Module : Data.Packed.Internal.Signatures
3-- Copyright : (c) Alberto Ruiz 2009
4-- License : BSD3
5-- Maintainer : Alberto Ruiz
6-- Stability : provisional
7--
8-- Signatures of the C functions.
9--
10
11
12module Data.Packed.Internal.Signatures where
13
14import Foreign.Ptr(Ptr)
15import Data.Complex(Complex)
16import Foreign.C.Types(CInt)
17
18type PF = Ptr Float --
19type PD = Ptr Double --
20type PQ = Ptr (Complex Float) --
21type PC = Ptr (Complex Double) --
22type TF = CInt -> PF -> IO CInt --
23type TFF = CInt -> PF -> TF --
24type TFV = CInt -> PF -> TV --
25type TVF = CInt -> PD -> TF --
26type TFFF = CInt -> PF -> TFF --
27type TV = CInt -> PD -> IO CInt --
28type TVV = CInt -> PD -> TV --
29type TVVV = CInt -> PD -> TVV --
30type TFM = CInt -> CInt -> PF -> IO CInt --
31type TFMFM = CInt -> CInt -> PF -> TFM --
32type TFMFMFM = CInt -> CInt -> PF -> TFMFM --
33type TM = CInt -> CInt -> PD -> IO CInt --
34type TMM = CInt -> CInt -> PD -> TM --
35type TVMM = CInt -> PD -> TMM --
36type TMVMM = CInt -> CInt -> PD -> TVMM --
37type TMMM = CInt -> CInt -> PD -> TMM --
38type TVM = CInt -> PD -> TM --
39type TVVM = CInt -> PD -> TVM --
40type TMV = CInt -> CInt -> PD -> TV --
41type TMMV = CInt -> CInt -> PD -> TMV --
42type TMVM = CInt -> CInt -> PD -> TVM --
43type TMMVM = CInt -> CInt -> PD -> TMVM --
44type TCM = CInt -> CInt -> PC -> IO CInt --
45type TCVCM = CInt -> PC -> TCM --
46type TCMCVCM = CInt -> CInt -> PC -> TCVCM --
47type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM --
48type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM --
49type TCMCM = CInt -> CInt -> PC -> TCM --
50type TVCM = CInt -> PD -> TCM --
51type TCMVCM = CInt -> CInt -> PC -> TVCM --
52type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM --
53type TCMCMCM = CInt -> CInt -> PC -> TCMCM --
54type TCV = CInt -> PC -> IO CInt --
55type TCVCV = CInt -> PC -> TCV --
56type TCVCVCV = CInt -> PC -> TCVCV --
57type TCVV = CInt -> PC -> TV --
58type TQV = CInt -> PQ -> IO CInt --
59type TQVQV = CInt -> PQ -> TQV --
60type TQVQVQV = CInt -> PQ -> TQVQV --
61type TQVF = CInt -> PQ -> TF --
62type TQM = CInt -> CInt -> PQ -> IO CInt --
63type TQMQM = CInt -> CInt -> PQ -> TQM --
64type TQMQMQM = CInt -> CInt -> PQ -> TQMQM --
65type TCMCV = CInt -> CInt -> PC -> TCV --
66type TVCV = CInt -> PD -> TCV --
67type TCVM = CInt -> PC -> TM --
68type TMCVM = CInt -> CInt -> PD -> TCVM --
69type TMMCVM = CInt -> CInt -> PD -> TMCVM --
70
diff --git a/packages/base/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs
deleted file mode 100644
index d0bc143..0000000
--- a/packages/base/src/Data/Packed/Internal/Vector.hs
+++ /dev/null
@@ -1,471 +0,0 @@
1{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-}
2-- |
3-- Module : Data.Packed.Internal.Vector
4-- Copyright : (c) Alberto Ruiz 2007
5-- License : BSD3
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-- Vector implementation
10--
11--------------------------------------------------------------------------------
12
13module Data.Packed.Internal.Vector (
14 Vector, dim,
15 fromList, toList, (|>),
16 vjoin, (@>), safe, at, at', subVector, takesV,
17 mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith,
18 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
19 foldVector, foldVectorG, foldLoop, foldVectorWithIndex,
20 createVector, vec,
21 asComplex, asReal, float2DoubleV, double2FloatV,
22 stepF, stepD, condF, condD,
23 conjugateQ, conjugateC,
24 cloneVector,
25 unsafeToForeignPtr,
26 unsafeFromForeignPtr,
27 unsafeWith
28) where
29
30import Data.Packed.Internal.Common
31import Data.Packed.Internal.Signatures
32import Foreign.Marshal.Array(peekArray, copyArray, advancePtr)
33import Foreign.ForeignPtr(ForeignPtr, castForeignPtr)
34import Foreign.Ptr(Ptr)
35import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf)
36import Foreign.C.Types
37import Data.Complex
38import Control.Monad(when)
39import System.IO.Unsafe(unsafePerformIO)
40
41#if __GLASGOW_HASKELL__ >= 605
42import GHC.ForeignPtr (mallocPlainForeignPtrBytes)
43#else
44import Foreign.ForeignPtr (mallocForeignPtrBytes)
45#endif
46
47import GHC.Base
48#if __GLASGOW_HASKELL__ < 612
49import GHC.IOBase hiding (liftIO)
50#endif
51
52import qualified Data.Vector.Storable as Vector
53import Data.Vector.Storable(Vector,
54 fromList,
55 unsafeToForeignPtr,
56 unsafeFromForeignPtr,
57 unsafeWith)
58
59
60-- | Number of elements
61dim :: (Storable t) => Vector t -> Int
62dim = Vector.length
63
64
65-- C-Haskell vector adapter
66-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
67vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
68vec x f = unsafeWith x $ \p -> do
69 let v g = do
70 g (fi $ dim x) p
71 f v
72{-# INLINE vec #-}
73
74
75-- allocates memory for a new vector
76createVector :: Storable a => Int -> IO (Vector a)
77createVector n = do
78 when (n < 0) $ error ("trying to createVector of negative dim: "++show n)
79 fp <- doMalloc undefined
80 return $ unsafeFromForeignPtr fp 0 n
81 where
82 --
83 -- Use the much cheaper Haskell heap allocated storage
84 -- for foreign pointer space we control
85 --
86 doMalloc :: Storable b => b -> IO (ForeignPtr b)
87 doMalloc dummy = do
88#if __GLASGOW_HASKELL__ >= 605
89 mallocPlainForeignPtrBytes (n * sizeOf dummy)
90#else
91 mallocForeignPtrBytes (n * sizeOf dummy)
92#endif
93
94{- | creates a Vector from a list:
95
96@> fromList [2,3,5,7]
974 |> [2.0,3.0,5.0,7.0]@
98
99-}
100
101safeRead v = inlinePerformIO . unsafeWith v
102{-# INLINE safeRead #-}
103
104inlinePerformIO :: IO a -> a
105inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
106{-# INLINE inlinePerformIO #-}
107
108{- | extracts the Vector elements to a list
109
110>>> toList (linspace 5 (1,10))
111[1.0,3.25,5.5,7.75,10.0]
112
113-}
114toList :: Storable a => Vector a -> [a]
115toList v = safeRead v $ peekArray (dim v)
116
117{- | Create a vector from a list of elements and explicit dimension. The input
118 list is explicitly truncated if it is too long, so it may safely
119 be used, for instance, with infinite lists.
120
121>>> 5 |> [1..]
122fromList [1.0,2.0,3.0,4.0,5.0]
123
124-}
125(|>) :: (Storable a) => Int -> [a] -> Vector a
126infixl 9 |>
127n |> l = if length l' == n
128 then fromList l'
129 else error "list too short for |>"
130 where l' = take n l
131
132
133-- | access to Vector elements without range checking
134at' :: Storable a => Vector a -> Int -> a
135at' v n = safeRead v $ flip peekElemOff n
136{-# INLINE at' #-}
137
138--
139-- turn off bounds checking with -funsafe at configure time.
140-- ghc will optimise away the salways true case at compile time.
141--
142#if defined(UNSAFE)
143safe :: Bool
144safe = False
145#else
146safe = True
147#endif
148
149-- | access to Vector elements with range checking.
150at :: Storable a => Vector a -> Int -> a
151at v n
152 | safe = if n >= 0 && n < dim v
153 then at' v n
154 else error "vector index out of range"
155 | otherwise = at' v n
156{-# INLINE at #-}
157
158{- | takes a number of consecutive elements from a Vector
159
160>>> subVector 2 3 (fromList [1..10])
161fromList [3.0,4.0,5.0]
162
163-}
164subVector :: Storable t => Int -- ^ index of the starting element
165 -> Int -- ^ number of elements to extract
166 -> Vector t -- ^ source
167 -> Vector t -- ^ result
168subVector = Vector.slice
169
170
171{- | Reads a vector position:
172
173>>> fromList [0..9] @> 7
1747.0
175
176-}
177(@>) :: Storable t => Vector t -> Int -> t
178infixl 9 @>
179(@>) = at
180
181
182{- | concatenate a list of vectors
183
184>>> vjoin [fromList [1..5::Double], konst 1 3]
185fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
186
187-}
188vjoin :: Storable t => [Vector t] -> Vector t
189vjoin [] = fromList []
190vjoin [v] = v
191vjoin as = unsafePerformIO $ do
192 let tot = sum (map dim as)
193 r <- createVector tot
194 unsafeWith r $ \ptr ->
195 joiner as tot ptr
196 return r
197 where joiner [] _ _ = return ()
198 joiner (v:cs) _ p = do
199 let n = dim v
200 unsafeWith v $ \pb -> copyArray p pb n
201 joiner cs 0 (advancePtr p n)
202
203
204{- | Extract consecutive subvectors of the given sizes.
205
206>>> takesV [3,4] (linspace 10 (1,10::Double))
207[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]
208
209-}
210takesV :: Storable t => [Int] -> Vector t -> [Vector t]
211takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w)
212 | otherwise = go ms w
213 where go [] _ = []
214 go (n:ns) v = subVector 0 n v
215 : go ns (subVector n (dim v - n) v)
216
217---------------------------------------------------------------
218
219-- | transforms a complex vector into a real vector with alternating real and imaginary parts
220asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a
221asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n)
222 where (fp,i,n) = unsafeToForeignPtr v
223
224-- | transforms a real vector into a complex vector with alternating real and imaginary parts
225asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a)
226asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2)
227 where (fp,i,n) = unsafeToForeignPtr v
228
229---------------------------------------------------------------
230
231float2DoubleV :: Vector Float -> Vector Double
232float2DoubleV v = unsafePerformIO $ do
233 r <- createVector (dim v)
234 app2 c_float2double vec v vec r "float2double"
235 return r
236
237double2FloatV :: Vector Double -> Vector Float
238double2FloatV v = unsafePerformIO $ do
239 r <- createVector (dim v)
240 app2 c_double2float vec v vec r "double2float2"
241 return r
242
243
244foreign import ccall unsafe "float2double" c_float2double:: TFV
245foreign import ccall unsafe "double2float" c_double2float:: TVF
246
247---------------------------------------------------------------
248
249stepF :: Vector Float -> Vector Float
250stepF v = unsafePerformIO $ do
251 r <- createVector (dim v)
252 app2 c_stepF vec v vec r "stepF"
253 return r
254
255stepD :: Vector Double -> Vector Double
256stepD v = unsafePerformIO $ do
257 r <- createVector (dim v)
258 app2 c_stepD vec v vec r "stepD"
259 return r
260
261foreign import ccall unsafe "stepF" c_stepF :: TFF
262foreign import ccall unsafe "stepD" c_stepD :: TVV
263
264---------------------------------------------------------------
265
266condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float
267condF x y l e g = unsafePerformIO $ do
268 r <- createVector (dim x)
269 app6 c_condF vec x vec y vec l vec e vec g vec r "condF"
270 return r
271
272condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double
273condD x y l e g = unsafePerformIO $ do
274 r <- createVector (dim x)
275 app6 c_condD vec x vec y vec l vec e vec g vec r "condD"
276 return r
277
278foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF
279foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV
280
281--------------------------------------------------------------------------------
282
283conjugateAux fun x = unsafePerformIO $ do
284 v <- createVector (dim x)
285 app2 fun vec x vec v "conjugateAux"
286 return v
287
288conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
289conjugateQ = conjugateAux c_conjugateQ
290foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV
291
292conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
293conjugateC = conjugateAux c_conjugateC
294foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV
295
296--------------------------------------------------------------------------------
297
298cloneVector :: Storable t => Vector t -> IO (Vector t)
299cloneVector v = do
300 let n = dim v
301 r <- createVector n
302 let f _ s _ d = copyArray d s n >> return 0
303 app2 f vec v vec r "cloneVector"
304 return r
305
306------------------------------------------------------------------
307
308-- | map on Vectors
309mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
310mapVector f v = unsafePerformIO $ do
311 w <- createVector (dim v)
312 unsafeWith v $ \p ->
313 unsafeWith w $ \q -> do
314 let go (-1) = return ()
315 go !k = do x <- peekElemOff p k
316 pokeElemOff q k (f x)
317 go (k-1)
318 go (dim v -1)
319 return w
320{-# INLINE mapVector #-}
321
322-- | zipWith for Vectors
323zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
324zipVectorWith f u v = unsafePerformIO $ do
325 let n = min (dim u) (dim v)
326 w <- createVector n
327 unsafeWith u $ \pu ->
328 unsafeWith v $ \pv ->
329 unsafeWith w $ \pw -> do
330 let go (-1) = return ()
331 go !k = do x <- peekElemOff pu k
332 y <- peekElemOff pv k
333 pokeElemOff pw k (f x y)
334 go (k-1)
335 go (n -1)
336 return w
337{-# INLINE zipVectorWith #-}
338
339-- | unzipWith for Vectors
340unzipVectorWith :: (Storable (a,b), Storable c, Storable d)
341 => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d)
342unzipVectorWith f u = unsafePerformIO $ do
343 let n = dim u
344 v <- createVector n
345 w <- createVector n
346 unsafeWith u $ \pu ->
347 unsafeWith v $ \pv ->
348 unsafeWith w $ \pw -> do
349 let go (-1) = return ()
350 go !k = do z <- peekElemOff pu k
351 let (x,y) = f z
352 pokeElemOff pv k x
353 pokeElemOff pw k y
354 go (k-1)
355 go (n-1)
356 return (v,w)
357{-# INLINE unzipVectorWith #-}
358
359foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b
360foldVector f x v = unsafePerformIO $
361 unsafeWith v $ \p -> do
362 let go (-1) s = return s
363 go !k !s = do y <- peekElemOff p k
364 go (k-1::Int) (f y s)
365 go (dim v -1) x
366{-# INLINE foldVector #-}
367
368-- the zero-indexed index is passed to the folding function
369foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b
370foldVectorWithIndex f x v = unsafePerformIO $
371 unsafeWith v $ \p -> do
372 let go (-1) s = return s
373 go !k !s = do y <- peekElemOff p k
374 go (k-1::Int) (f k y s)
375 go (dim v -1) x
376{-# INLINE foldVectorWithIndex #-}
377
378foldLoop f s0 d = go (d - 1) s0
379 where
380 go 0 s = f (0::Int) s
381 go !j !s = go (j - 1) (f j s)
382
383foldVectorG f s0 v = foldLoop g s0 (dim v)
384 where g !k !s = f k (at' v) s
385 {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
386{-# INLINE foldVectorG #-}
387
388-------------------------------------------------------------------
389
390-- | monadic map over Vectors
391-- the monad @m@ must be strict
392mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b)
393mapVectorM f v = do
394 w <- return $! unsafePerformIO $! createVector (dim v)
395 mapVectorM' w 0 (dim v -1)
396 return w
397 where mapVectorM' w' !k !t
398 | k == t = do
399 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
400 y <- f x
401 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
402 | otherwise = do
403 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
404 y <- f x
405 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
406 mapVectorM' w' (k+1) t
407{-# INLINE mapVectorM #-}
408
409-- | monadic map over Vectors
410mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m ()
411mapVectorM_ f v = do
412 mapVectorM' 0 (dim v -1)
413 where mapVectorM' !k !t
414 | k == t = do
415 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
416 f x
417 | otherwise = do
418 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
419 _ <- f x
420 mapVectorM' (k+1) t
421{-# INLINE mapVectorM_ #-}
422
423-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
424-- the monad @m@ must be strict
425mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b)
426mapVectorWithIndexM f v = do
427 w <- return $! unsafePerformIO $! createVector (dim v)
428 mapVectorM' w 0 (dim v -1)
429 return w
430 where mapVectorM' w' !k !t
431 | k == t = do
432 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
433 y <- f k x
434 return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
435 | otherwise = do
436 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
437 y <- f k x
438 _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y
439 mapVectorM' w' (k+1) t
440{-# INLINE mapVectorWithIndexM #-}
441
442-- | monadic map over Vectors with the zero-indexed index passed to the mapping function
443mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m ()
444mapVectorWithIndexM_ f v = do
445 mapVectorM' 0 (dim v -1)
446 where mapVectorM' !k !t
447 | k == t = do
448 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
449 f k x
450 | otherwise = do
451 x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k
452 _ <- f k x
453 mapVectorM' (k+1) t
454{-# INLINE mapVectorWithIndexM_ #-}
455
456
457mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b
458--mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b])
459mapVectorWithIndex f v = unsafePerformIO $ do
460 w <- createVector (dim v)
461 unsafeWith v $ \p ->
462 unsafeWith w $ \q -> do
463 let go (-1) = return ()
464 go !k = do x <- peekElemOff p k
465 pokeElemOff q k (f k x)
466 go (k-1)
467 go (dim v -1)
468 return w
469{-# INLINE mapVectorWithIndex #-}
470
471
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
deleted file mode 100644
index 70b9232..0000000
--- a/packages/base/src/Data/Packed/Matrix.hs
+++ /dev/null
@@ -1,494 +0,0 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-}
6{-# LANGUAGE CPP #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Data.Packed.Matrix
11-- Copyright : (c) Alberto Ruiz 2007-10
12-- License : BSD3
13-- Maintainer : Alberto Ruiz
14-- Stability : provisional
15--
16-- A Matrix representation suitable for numerical computations using LAPACK and GSL.
17--
18-- This module provides basic functions for manipulation of structure.
19
20-----------------------------------------------------------------------------
21{-# OPTIONS_HADDOCK hide #-}
22
23module Data.Packed.Matrix (
24 Matrix,
25 Element,
26 rows,cols,
27 (><),
28 trans,
29 reshape, flatten,
30 fromLists, toLists, buildMatrix,
31 (@@>),
32 asRow, asColumn,
33 fromRows, toRows, fromColumns, toColumns,
34 fromBlocks, diagBlock, toBlocks, toBlocksEvery,
35 repmat,
36 flipud, fliprl,
37 subMatrix, takeRows, dropRows, takeColumns, dropColumns,
38 extractRows, extractColumns,
39 diagRect, takeDiag,
40 mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
41 liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D
42) where
43
44import Data.Packed.Internal
45import qualified Data.Packed.ST as ST
46import Data.Array
47
48import Data.List(transpose,intersperse)
49import Foreign.Storable(Storable)
50import Control.Monad(liftM)
51
52-------------------------------------------------------------------
53
54#ifdef BINARY
55
56import Data.Binary
57
58instance (Binary (Vector a), Element a) => Binary (Matrix a) where
59 put m = do
60 put (cols m)
61 put (flatten m)
62 get = do
63 c <- get
64 v <- get
65 return (reshape c v)
66
67#endif
68
69-------------------------------------------------------------------
70
71instance (Show a, Element a) => (Show (Matrix a)) where
72 show m | rows m == 0 || cols m == 0 = sizes m ++" []"
73 show m = (sizes m++) . dsp . map (map show) . toLists $ m
74
75sizes m = "("++show (rows m)++"><"++show (cols m)++")\n"
76
77dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
78 where
79 mt = transpose as
80 longs = map (maximum . map length) mt
81 mtp = zipWith (\a b -> map (pad a) b) longs mt
82 pad n str = replicate (n - length str) ' ' ++ str
83 unwords' = concat . intersperse ", "
84
85------------------------------------------------------------------
86
87instance (Element a, Read a) => Read (Matrix a) where
88 readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
89 where (thing,rest) = breakAt ']' s
90 (dims,listnums) = breakAt ')' thing
91 cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
92 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
93
94
95breakAt c l = (a++[c],tail b) where
96 (a,b) = break (==c) l
97
98------------------------------------------------------------------
99
100-- | creates a matrix from a vertical list of matrices
101joinVert :: Element t => [Matrix t] -> Matrix t
102joinVert [] = emptyM 0 0
103joinVert ms = case common cols ms of
104 Nothing -> error "(impossible) joinVert on matrices with different number of columns"
105 Just c -> matrixFromVector RowMajor (sum (map rows ms)) c $ vjoin (map flatten ms)
106
107-- | creates a matrix from a horizontal list of matrices
108joinHoriz :: Element t => [Matrix t] -> Matrix t
109joinHoriz ms = trans. joinVert . map trans $ ms
110
111{- | Create a matrix from blocks given as a list of lists of matrices.
112
113Single row-column components are automatically expanded to match the
114corresponding common row and column:
115
116@
117disp = putStr . dispf 2
118@
119
120>>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]]
1218x10
1221 0 0 0 0 7 7 7 10 20
1230 1 0 0 0 7 7 7 10 20
1240 0 1 0 0 7 7 7 10 20
1250 0 0 1 0 7 7 7 10 20
1260 0 0 0 1 7 7 7 10 20
1273 3 3 3 3 1 0 0 0 0
1283 3 3 3 3 0 2 0 0 0
1293 3 3 3 3 0 0 3 0 0
130
131-}
132fromBlocks :: Element t => [[Matrix t]] -> Matrix t
133fromBlocks = fromBlocksRaw . adaptBlocks
134
135fromBlocksRaw mms = joinVert . map joinHoriz $ mms
136
137adaptBlocks ms = ms' where
138 bc = case common length ms of
139 Just c -> c
140 Nothing -> error "fromBlocks requires rectangular [[Matrix]]"
141 rs = map (compatdim . map rows) ms
142 cs = map (compatdim . map cols) (transpose ms)
143 szs = sequence [rs,cs]
144 ms' = splitEvery bc $ zipWith g szs (concat ms)
145
146 g [Just nr,Just nc] m
147 | nr == r && nc == c = m
148 | r == 1 && c == 1 = matrixFromVector RowMajor nr nc (constantD x (nr*nc))
149 | r == 1 = fromRows (replicate nr (flatten m))
150 | otherwise = fromColumns (replicate nc (flatten m))
151 where
152 r = rows m
153 c = cols m
154 x = m@@>(0,0)
155 g _ _ = error "inconsistent dimensions in fromBlocks"
156
157
158--------------------------------------------------------------------------------
159
160{- | create a block diagonal matrix
161
162>>> disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]]
1637x8
1641 1 0 0 0 0 0 0
1651 1 0 0 0 0 0 0
1660 0 2 2 2 2 2 0
1670 0 2 2 2 2 2 0
1680 0 2 2 2 2 2 0
1690 0 0 0 0 0 0 5
1700 0 0 0 0 0 0 7
171
172>>> diagBlock [(0><4)[], konst 2 (2,3)] :: Matrix Double
173(2><7)
174 [ 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0
175 , 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 ]
176
177-}
178diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
179diagBlock ms = fromBlocks $ zipWith f ms [0..]
180 where
181 f m k = take n $ replicate k z ++ m : repeat z
182 n = length ms
183 z = (1><1) [0]
184
185--------------------------------------------------------------------------------
186
187
188-- | Reverse rows
189flipud :: Element t => Matrix t -> Matrix t
190flipud m = extractRows [r-1,r-2 .. 0] $ m
191 where
192 r = rows m
193
194-- | Reverse columns
195fliprl :: Element t => Matrix t -> Matrix t
196fliprl m = extractColumns [c-1,c-2 .. 0] $ m
197 where
198 c = cols m
199
200------------------------------------------------------------
201
202{- | creates a rectangular diagonal matrix:
203
204>>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
205(4><5)
206 [ 10.0, 7.0, 7.0, 7.0, 7.0
207 , 7.0, 20.0, 7.0, 7.0, 7.0
208 , 7.0, 7.0, 30.0, 7.0, 7.0
209 , 7.0, 7.0, 7.0, 7.0, 7.0 ]
210
211-}
212diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
213diagRect z v r c = ST.runSTMatrix $ do
214 m <- ST.newMatrix z r c
215 let d = min r c `min` (dim v)
216 mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
217 return m
218
219-- | extracts the diagonal from a rectangular matrix
220takeDiag :: (Element t) => Matrix t -> Vector t
221takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]]
222
223------------------------------------------------------------
224
225{- | create a general matrix
226
227>>> (2><3) [2, 4, 7+2*𝑖, -3, 11, 0]
228(2><3)
229 [ 2.0 :+ 0.0, 4.0 :+ 0.0, 7.0 :+ 2.0
230 , (-3.0) :+ (-0.0), 11.0 :+ 0.0, 0.0 :+ 0.0 ]
231
232The input list is explicitly truncated, so that it can
233safely be used with lists that are too long (like infinite lists).
234
235>>> (2><3)[1..]
236(2><3)
237 [ 1.0, 2.0, 3.0
238 , 4.0, 5.0, 6.0 ]
239
240This is the format produced by the instances of Show (Matrix a), which
241can also be used for input.
242
243-}
244(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
245r >< c = f where
246 f l | dim v == r*c = matrixFromVector RowMajor r c v
247 | otherwise = error $ "inconsistent list size = "
248 ++show (dim v) ++" in ("++show r++"><"++show c++")"
249 where v = fromList $ take (r*c) l
250
251----------------------------------------------------------------
252
253-- | Creates a matrix with the first n rows of another matrix
254takeRows :: Element t => Int -> Matrix t -> Matrix t
255takeRows n mt = subMatrix (0,0) (n, cols mt) mt
256-- | Creates a copy of a matrix without the first n rows
257dropRows :: Element t => Int -> Matrix t -> Matrix t
258dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
259-- |Creates a matrix with the first n columns of another matrix
260takeColumns :: Element t => Int -> Matrix t -> Matrix t
261takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
262-- | Creates a copy of a matrix without the first n columns
263dropColumns :: Element t => Int -> Matrix t -> Matrix t
264dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
265
266----------------------------------------------------------------
267
268{- | Creates a 'Matrix' from a list of lists (considered as rows).
269
270>>> fromLists [[1,2],[3,4],[5,6]]
271(3><2)
272 [ 1.0, 2.0
273 , 3.0, 4.0
274 , 5.0, 6.0 ]
275
276-}
277fromLists :: Element t => [[t]] -> Matrix t
278fromLists = fromRows . map fromList
279
280-- | creates a 1-row matrix from a vector
281--
282-- >>> asRow (fromList [1..5])
283-- (1><5)
284-- [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
285--
286asRow :: Storable a => Vector a -> Matrix a
287asRow = trans . asColumn
288
289-- | creates a 1-column matrix from a vector
290--
291-- >>> asColumn (fromList [1..5])
292-- (5><1)
293-- [ 1.0
294-- , 2.0
295-- , 3.0
296-- , 4.0
297-- , 5.0 ]
298--
299asColumn :: Storable a => Vector a -> Matrix a
300asColumn v = reshape 1 v
301
302
303
304{- | creates a Matrix of the specified size using the supplied function to
305 to map the row\/column position to the value at that row\/column position.
306
307@> buildMatrix 3 4 (\\(r,c) -> fromIntegral r * fromIntegral c)
308(3><4)
309 [ 0.0, 0.0, 0.0, 0.0, 0.0
310 , 0.0, 1.0, 2.0, 3.0, 4.0
311 , 0.0, 2.0, 4.0, 6.0, 8.0]@
312
313Hilbert matrix of order N:
314
315@hilb n = buildMatrix n n (\\(i,j)->1/(fromIntegral i + fromIntegral j +1))@
316
317-}
318buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a
319buildMatrix rc cc f =
320 fromLists $ map (map f)
321 $ map (\ ri -> map (\ ci -> (ri, ci)) [0 .. (cc - 1)]) [0 .. (rc - 1)]
322
323-----------------------------------------------------
324
325fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e
326fromArray2D m = (r><c) (elems m)
327 where ((r0,c0),(r1,c1)) = bounds m
328 r = r1-r0+1
329 c = c1-c0+1
330
331
332-- | rearranges the rows of a matrix according to the order given in a list of integers.
333extractRows :: Element t => [Int] -> Matrix t -> Matrix t
334extractRows [] m = emptyM 0 (cols m)
335extractRows l m = fromRows $ extract (toRows m) l
336 where
337 extract l' is = [l'!!i | i<- map verify is]
338 verify k
339 | k >= 0 && k < rows m = k
340 | otherwise = error $ "can't extract row "
341 ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
342
343-- | rearranges the rows of a matrix according to the order given in a list of integers.
344extractColumns :: Element t => [Int] -> Matrix t -> Matrix t
345extractColumns l m = trans . extractRows (map verify l) . trans $ m
346 where
347 verify k
348 | k >= 0 && k < cols m = k
349 | otherwise = error $ "can't extract column "
350 ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
351
352
353
354{- | creates matrix by repetition of a matrix a given number of rows and columns
355
356>>> repmat (ident 2) 2 3
357(4><6)
358 [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
359 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
360 , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
361 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]
362
363-}
364repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
365repmat m r c
366 | r == 0 || c == 0 = emptyM (r*rows m) (c*cols m)
367 | otherwise = fromBlocks $ replicate r $ replicate c $ m
368
369-- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.
370liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
371liftMatrix2Auto f m1 m2
372 | compat' m1 m2 = lM f m1 m2
373 | ok = lM f m1' m2'
374 | otherwise = error $ "nonconformable matrices in liftMatrix2Auto: " ++ shSize m1 ++ ", " ++ shSize m2
375 where
376 (r1,c1) = size m1
377 (r2,c2) = size m2
378 r = max r1 r2
379 c = max c1 c2
380 r0 = min r1 r2
381 c0 = min c1 c2
382 ok = r0 == 1 || r1 == r2 && c0 == 1 || c1 == c2
383 m1' = conformMTo (r,c) m1
384 m2' = conformMTo (r,c) m2
385
386-- FIXME do not flatten if equal order
387lM f m1 m2 = matrixFromVector
388 RowMajor
389 (max (rows m1) (rows m2))
390 (max (cols m1) (cols m2))
391 (f (flatten m1) (flatten m2))
392
393compat' :: Matrix a -> Matrix b -> Bool
394compat' m1 m2 = s1 == (1,1) || s2 == (1,1) || s1 == s2
395 where
396 s1 = size m1
397 s2 = size m2
398
399------------------------------------------------------------
400
401toBlockRows [r] m
402 | r == rows m = [m]
403toBlockRows rs m
404 | cols m > 0 = map (reshape (cols m)) (takesV szs (flatten m))
405 | otherwise = map g rs
406 where
407 szs = map (* cols m) rs
408 g k = (k><0)[]
409
410toBlockCols [c] m | c == cols m = [m]
411toBlockCols cs m = map trans . toBlockRows cs . trans $ m
412
413-- | Partition a matrix into blocks with the given numbers of rows and columns.
414-- The remaining rows and columns are discarded.
415toBlocks :: (Element t) => [Int] -> [Int] -> Matrix t -> [[Matrix t]]
416toBlocks rs cs m
417 | ok = map (toBlockCols cs) . toBlockRows rs $ m
418 | otherwise = error $ "toBlocks: bad partition: "++show rs++" "++show cs
419 ++ " "++shSize m
420 where
421 ok = sum rs <= rows m && sum cs <= cols m && all (>=0) rs && all (>=0) cs
422
423-- | Fully partition a matrix into blocks of the same size. If the dimensions are not
424-- a multiple of the given size the last blocks will be smaller.
425toBlocksEvery :: (Element t) => Int -> Int -> Matrix t -> [[Matrix t]]
426toBlocksEvery r c m
427 | r < 1 || c < 1 = error $ "toBlocksEvery expects block sizes > 0, given "++show r++" and "++ show c
428 | otherwise = toBlocks rs cs m
429 where
430 (qr,rr) = rows m `divMod` r
431 (qc,rc) = cols m `divMod` c
432 rs = replicate qr r ++ if rr > 0 then [rr] else []
433 cs = replicate qc c ++ if rc > 0 then [rc] else []
434
435-------------------------------------------------------------------
436
437-- Given a column number and a function taking matrix indexes, returns
438-- a function which takes vector indexes (that can be used on the
439-- flattened matrix).
440mk :: Int -> ((Int, Int) -> t) -> (Int -> t)
441mk c g = \k -> g (divMod k c)
442
443{- |
444
445>>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..])
446m[0,0] = 1
447m[0,1] = 2
448m[0,2] = 3
449m[1,0] = 4
450m[1,1] = 5
451m[1,2] = 6
452
453-}
454mapMatrixWithIndexM_
455 :: (Element a, Num a, Monad m) =>
456 ((Int, Int) -> a -> m ()) -> Matrix a -> m ()
457mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
458 where
459 c = cols m
460
461{- |
462
463>>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
464Just (3><3)
465 [ 100.0, 1.0, 2.0
466 , 10.0, 111.0, 12.0
467 , 20.0, 21.0, 122.0 ]
468
469-}
470mapMatrixWithIndexM
471 :: (Element a, Storable b, Monad m) =>
472 ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
473mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m
474 where
475 c = cols m
476
477{- |
478
479>>> mapMatrixWithIndex (\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
480(3><3)
481 [ 100.0, 1.0, 2.0
482 , 10.0, 111.0, 12.0
483 , 20.0, 21.0, 122.0 ]
484
485 -}
486mapMatrixWithIndex
487 :: (Element a, Storable b) =>
488 ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
489mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
490 where
491 c = cols m
492
493mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b
494mapMatrix f = liftMatrix (mapVector f)
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
deleted file mode 100644
index 6027f43..0000000
--- a/packages/base/src/Data/Packed/Numeric.hs
+++ /dev/null
@@ -1,299 +0,0 @@
1{-# LANGUAGE FlexibleContexts #-}
2{-# LANGUAGE FlexibleInstances #-}
3{-# LANGUAGE MultiParamTypeClasses #-}
4{-# LANGUAGE FunctionalDependencies #-}
5{-# LANGUAGE UndecidableInstances #-}
6
7-----------------------------------------------------------------------------
8-- |
9-- Module : Data.Packed.Numeric
10-- Copyright : (c) Alberto Ruiz 2010-14
11-- License : BSD3
12-- Maintainer : Alberto Ruiz
13-- Stability : provisional
14--
15-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
16--
17-- The 'Container' class is used to define optimized generic functions which work
18-- on 'Vector' and 'Matrix' with real or complex elements.
19--
20-- Some of these functions are also available in the instances of the standard
21-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
22--
23-----------------------------------------------------------------------------
24{-# OPTIONS_HADDOCK hide #-}
25
26module Data.Packed.Numeric (
27 -- * Basic functions
28 module Data.Packed,
29 Konst(..), Build(..),
30 linspace,
31 diag, ident,
32 ctrans,
33 -- * Generic operations
34 Container(..), Numeric,
35 -- add, mul, sub, divide, equal, scaleRecip, addConstant,
36 scalar, conj, scale, arctan2, cmap,
37 atIndex, minIndex, maxIndex, minElement, maxElement,
38 sumElements, prodElements,
39 step, cond, find, assoc, accum,
40 Transposable(..), Linear(..),
41 -- * Matrix product
42 Product(..), udot, dot, (<·>), (#>), app,
43 Mul(..),
44 (<.>),
45 optimiseMult,
46 mXm,mXv,vXm,LSDiv,(<\>),
47 outer, kronecker,
48 -- * Random numbers
49 RandDist(..),
50 randomVector,
51 gaussianSample,
52 uniformSample,
53 meanCov,
54 -- * sorting
55 sortVector,
56 -- * Element conversion
57 Convert(..),
58 Complexable(),
59 RealElement(),
60 RealOf, ComplexOf, SingleOf, DoubleOf,
61 roundVector,
62 IndexOf,
63 module Data.Complex,
64 -- * IO
65 module Data.Packed.IO,
66 -- * Misc
67 Testable(..)
68) where
69
70import Data.Packed
71import Data.Packed.Internal.Numeric
72import Data.Complex
73import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
74import Data.Monoid(Monoid(mconcat))
75import Data.Packed.IO
76import Numeric.LinearAlgebra.Random
77
78------------------------------------------------------------------
79
80{- | Creates a real vector containing a range of values:
81
82>>> linspace 5 (-3,7::Double)
83fromList [-3.0,-0.5,2.0,4.5,7.0]@
84
85>>> linspace 5 (8,2+i) :: Vector (Complex Double)
86fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
87
88Logarithmic spacing can be defined as follows:
89
90@logspace n (a,b) = 10 ** linspace n (a,b)@
91-}
92linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
93linspace 0 _ = fromList[]
94linspace 1 (a,b) = fromList[(a+b)/2]
95linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
96 where s = (b-a)/fromIntegral (n-1)
97
98--------------------------------------------------------------------------------
99
100infixl 7 <.>
101-- | An infix synonym for 'dot'
102(<.>) :: Numeric t => Vector t -> Vector t -> t
103(<.>) = dot
104
105
106infixr 8 <·>, #>
107
108{- | infix synonym for 'dot'
109
110>>> vector [1,2,3,4] <·> vector [-2,0,1,1]
1115.0
112
113>>> let 𝑖 = 0:+1 :: ℂ
114>>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖]
1152.0 :+ 0.0
116
117(the dot symbol "·" is obtained by Alt-Gr .)
118
119-}
120(<·>) :: Numeric t => Vector t -> Vector t -> t
121(<·>) = dot
122
123
124{- | infix synonym for 'app'
125
126>>> let m = (2><3) [1..]
127>>> m
128(2><3)
129 [ 1.0, 2.0, 3.0
130 , 4.0, 5.0, 6.0 ]
131
132>>> let v = vector [10,20,30]
133
134>>> m #> v
135fromList [140.0,320.0]
136
137-}
138(#>) :: Numeric t => Matrix t -> Vector t -> Vector t
139(#>) = mXv
140
141-- | dense matrix-vector product
142app :: Numeric t => Matrix t -> Vector t -> Vector t
143app = (#>)
144
145--------------------------------------------------------------------------------
146
147class Mul a b c | a b -> c where
148 infixl 7 <>
149 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
150 (<>) :: Product t => a t -> b t -> c t
151
152instance Mul Matrix Matrix Matrix where
153 (<>) = mXm
154
155instance Mul Matrix Vector Vector where
156 (<>) m v = flatten $ m <> asColumn v
157
158instance Mul Vector Matrix Vector where
159 (<>) v m = flatten $ asRow v <> m
160
161--------------------------------------------------------------------------------
162
163{- | Least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
164
165@
166a = (3><2)
167 [ 1.0, 2.0
168 , 2.0, 4.0
169 , 2.0, -1.0 ]
170@
171
172@
173v = vector [13.0,27.0,1.0]
174@
175
176>>> let x = a <\> v
177>>> x
178fromList [3.0799999999999996,5.159999999999999]
179
180>>> a #> x
181fromList [13.399999999999999,26.799999999999997,1.0]
182
183It also admits multiple right-hand sides stored as columns in a matrix.
184
185-}
186infixl 7 <\>
187(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
188(<\>) = linSolve
189
190class LSDiv c
191 where
192 linSolve :: Field t => Matrix t -> c t -> c t
193
194instance LSDiv Vector
195 where
196 linSolve m v = flatten (linearSolveSVD m (reshape 1 v))
197
198instance LSDiv Matrix
199 where
200 linSolve = linearSolveSVD
201
202--------------------------------------------------------------------------------
203
204class Konst e d c | d -> c, c -> d
205 where
206 -- |
207 -- >>> konst 7 3 :: Vector Float
208 -- fromList [7.0,7.0,7.0]
209 --
210 -- >>> konst i (3::Int,4::Int)
211 -- (3><4)
212 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
213 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
214 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
215 --
216 konst :: e -> d -> c e
217
218instance Container Vector e => Konst e Int Vector
219 where
220 konst = konst'
221
222instance Container Vector e => Konst e (Int,Int) Matrix
223 where
224 konst = konst'
225
226--------------------------------------------------------------------------------
227
228class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
229 where
230 -- |
231 -- >>> build 5 (**2) :: Vector Double
232 -- fromList [0.0,1.0,4.0,9.0,16.0]
233 --
234 -- Hilbert matrix of order N:
235 --
236 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
237 -- >>> putStr . dispf 2 $ hilb 3
238 -- 3x3
239 -- 1.00 0.50 0.33
240 -- 0.50 0.33 0.25
241 -- 0.33 0.25 0.20
242 --
243 build :: d -> f -> c e
244
245instance Container Vector e => Build Int (e -> e) Vector e
246 where
247 build = build'
248
249instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
250 where
251 build = build'
252
253--------------------------------------------------------------------------------
254
255-- @dot u v = 'udot' ('conj' u) v@
256dot :: (Numeric t) => Vector t -> Vector t -> t
257dot u v = udot (conj u) v
258
259--------------------------------------------------------------------------------
260
261optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
262optimiseMult = mconcat
263
264--------------------------------------------------------------------------------
265
266
267{- | Compute mean vector and covariance matrix of the rows of a matrix.
268
269>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
270(fromList [4.010341078059521,5.0197204699640405],
271(2><2)
272 [ 1.9862461923890056, -1.0127225830525157e-2
273 , -1.0127225830525157e-2, 3.0373954915729318 ])
274
275-}
276meanCov :: Matrix Double -> (Vector Double, Matrix Double)
277meanCov x = (med,cov) where
278 r = rows x
279 k = 1 / fromIntegral r
280 med = konst k r `vXm` x
281 meds = konst 1 r `outer` med
282 xc = x `sub` meds
283 cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)
284
285--------------------------------------------------------------------------------
286
287class ( Container Vector t
288 , Container Matrix t
289 , Konst t Int Vector
290 , Konst t (Int,Int) Matrix
291 , Product t
292 ) => Numeric t
293
294instance Numeric Double
295instance Numeric (Complex Double)
296instance Numeric Float
297instance Numeric (Complex Float)
298
299
diff --git a/packages/base/src/Data/Packed/ST.hs b/packages/base/src/Data/Packed/ST.hs
deleted file mode 100644
index 5c45c7b..0000000
--- a/packages/base/src/Data/Packed/ST.hs
+++ /dev/null
@@ -1,178 +0,0 @@
1{-# LANGUAGE CPP #-}
2{-# LANGUAGE TypeOperators #-}
3{-# LANGUAGE Rank2Types #-}
4{-# LANGUAGE BangPatterns #-}
5-----------------------------------------------------------------------------
6-- |
7-- Module : Data.Packed.ST
8-- Copyright : (c) Alberto Ruiz 2008
9-- License : BSD3
10-- Maintainer : Alberto Ruiz
11-- Stability : provisional
12--
13-- In-place manipulation inside the ST monad.
14-- See examples/inplace.hs in the distribution.
15--
16-----------------------------------------------------------------------------
17{-# OPTIONS_HADDOCK hide #-}
18
19module Data.Packed.ST (
20 -- * Mutable Vectors
21 STVector, newVector, thawVector, freezeVector, runSTVector,
22 readVector, writeVector, modifyVector, liftSTVector,
23 -- * Mutable Matrices
24 STMatrix, newMatrix, thawMatrix, freezeMatrix, runSTMatrix,
25 readMatrix, writeMatrix, modifyMatrix, liftSTMatrix,
26 -- * Unsafe functions
27 newUndefinedVector,
28 unsafeReadVector, unsafeWriteVector,
29 unsafeThawVector, unsafeFreezeVector,
30 newUndefinedMatrix,
31 unsafeReadMatrix, unsafeWriteMatrix,
32 unsafeThawMatrix, unsafeFreezeMatrix
33) where
34
35import Data.Packed.Internal
36
37import Control.Monad.ST(ST, runST)
38import Foreign.Storable(Storable, peekElemOff, pokeElemOff)
39
40#if MIN_VERSION_base(4,4,0)
41import Control.Monad.ST.Unsafe(unsafeIOToST)
42#else
43import Control.Monad.ST(unsafeIOToST)
44#endif
45
46{-# INLINE ioReadV #-}
47ioReadV :: Storable t => Vector t -> Int -> IO t
48ioReadV v k = unsafeWith v $ \s -> peekElemOff s k
49
50{-# INLINE ioWriteV #-}
51ioWriteV :: Storable t => Vector t -> Int -> t -> IO ()
52ioWriteV v k x = unsafeWith v $ \s -> pokeElemOff s k x
53
54newtype STVector s t = STVector (Vector t)
55
56thawVector :: Storable t => Vector t -> ST s (STVector s t)
57thawVector = unsafeIOToST . fmap STVector . cloneVector
58
59unsafeThawVector :: Storable t => Vector t -> ST s (STVector s t)
60unsafeThawVector = unsafeIOToST . return . STVector
61
62runSTVector :: Storable t => (forall s . ST s (STVector s t)) -> Vector t
63runSTVector st = runST (st >>= unsafeFreezeVector)
64
65{-# INLINE unsafeReadVector #-}
66unsafeReadVector :: Storable t => STVector s t -> Int -> ST s t
67unsafeReadVector (STVector x) = unsafeIOToST . ioReadV x
68
69{-# INLINE unsafeWriteVector #-}
70unsafeWriteVector :: Storable t => STVector s t -> Int -> t -> ST s ()
71unsafeWriteVector (STVector x) k = unsafeIOToST . ioWriteV x k
72
73{-# INLINE modifyVector #-}
74modifyVector :: (Storable t) => STVector s t -> Int -> (t -> t) -> ST s ()
75modifyVector x k f = readVector x k >>= return . f >>= unsafeWriteVector x k
76
77liftSTVector :: (Storable t) => (Vector t -> a) -> STVector s1 t -> ST s2 a
78liftSTVector f (STVector x) = unsafeIOToST . fmap f . cloneVector $ x
79
80freezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
81freezeVector v = liftSTVector id v
82
83unsafeFreezeVector :: (Storable t) => STVector s1 t -> ST s2 (Vector t)
84unsafeFreezeVector (STVector x) = unsafeIOToST . return $ x
85
86{-# INLINE safeIndexV #-}
87safeIndexV f (STVector v) k
88 | k < 0 || k>= dim v = error $ "out of range error in vector (dim="
89 ++show (dim v)++", pos="++show k++")"
90 | otherwise = f (STVector v) k
91
92{-# INLINE readVector #-}
93readVector :: Storable t => STVector s t -> Int -> ST s t
94readVector = safeIndexV unsafeReadVector
95
96{-# INLINE writeVector #-}
97writeVector :: Storable t => STVector s t -> Int -> t -> ST s ()
98writeVector = safeIndexV unsafeWriteVector
99
100newUndefinedVector :: Storable t => Int -> ST s (STVector s t)
101newUndefinedVector = unsafeIOToST . fmap STVector . createVector
102
103{-# INLINE newVector #-}
104newVector :: Storable t => t -> Int -> ST s (STVector s t)
105newVector x n = do
106 v <- newUndefinedVector n
107 let go (-1) = return v
108 go !k = unsafeWriteVector v k x >> go (k-1 :: Int)
109 go (n-1)
110
111-------------------------------------------------------------------------
112
113{-# INLINE ioReadM #-}
114ioReadM :: Storable t => Matrix t -> Int -> Int -> IO t
115ioReadM (Matrix _ nc cv RowMajor) r c = ioReadV cv (r*nc+c)
116ioReadM (Matrix nr _ fv ColumnMajor) r c = ioReadV fv (c*nr+r)
117
118{-# INLINE ioWriteM #-}
119ioWriteM :: Storable t => Matrix t -> Int -> Int -> t -> IO ()
120ioWriteM (Matrix _ nc cv RowMajor) r c val = ioWriteV cv (r*nc+c) val
121ioWriteM (Matrix nr _ fv ColumnMajor) r c val = ioWriteV fv (c*nr+r) val
122
123newtype STMatrix s t = STMatrix (Matrix t)
124
125thawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
126thawMatrix = unsafeIOToST . fmap STMatrix . cloneMatrix
127
128unsafeThawMatrix :: Storable t => Matrix t -> ST s (STMatrix s t)
129unsafeThawMatrix = unsafeIOToST . return . STMatrix
130
131runSTMatrix :: Storable t => (forall s . ST s (STMatrix s t)) -> Matrix t
132runSTMatrix st = runST (st >>= unsafeFreezeMatrix)
133
134{-# INLINE unsafeReadMatrix #-}
135unsafeReadMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
136unsafeReadMatrix (STMatrix x) r = unsafeIOToST . ioReadM x r
137
138{-# INLINE unsafeWriteMatrix #-}
139unsafeWriteMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
140unsafeWriteMatrix (STMatrix x) r c = unsafeIOToST . ioWriteM x r c
141
142{-# INLINE modifyMatrix #-}
143modifyMatrix :: (Storable t) => STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
144modifyMatrix x r c f = readMatrix x r c >>= return . f >>= unsafeWriteMatrix x r c
145
146liftSTMatrix :: (Storable t) => (Matrix t -> a) -> STMatrix s1 t -> ST s2 a
147liftSTMatrix f (STMatrix x) = unsafeIOToST . fmap f . cloneMatrix $ x
148
149unsafeFreezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
150unsafeFreezeMatrix (STMatrix x) = unsafeIOToST . return $ x
151
152freezeMatrix :: (Storable t) => STMatrix s1 t -> ST s2 (Matrix t)
153freezeMatrix m = liftSTMatrix id m
154
155cloneMatrix (Matrix r c d o) = cloneVector d >>= return . (\d' -> Matrix r c d' o)
156
157{-# INLINE safeIndexM #-}
158safeIndexM f (STMatrix m) r c
159 | r<0 || r>=rows m ||
160 c<0 || c>=cols m = error $ "out of range error in matrix (size="
161 ++show (rows m,cols m)++", pos="++show (r,c)++")"
162 | otherwise = f (STMatrix m) r c
163
164{-# INLINE readMatrix #-}
165readMatrix :: Storable t => STMatrix s t -> Int -> Int -> ST s t
166readMatrix = safeIndexM unsafeReadMatrix
167
168{-# INLINE writeMatrix #-}
169writeMatrix :: Storable t => STMatrix s t -> Int -> Int -> t -> ST s ()
170writeMatrix = safeIndexM unsafeWriteMatrix
171
172newUndefinedMatrix :: Storable t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t)
173newUndefinedMatrix ord r c = unsafeIOToST $ fmap STMatrix $ createMatrix ord r c
174
175{-# NOINLINE newMatrix #-}
176newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t)
177newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c)
178
diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs
deleted file mode 100644
index 2104f52..0000000
--- a/packages/base/src/Data/Packed/Vector.hs
+++ /dev/null
@@ -1,125 +0,0 @@
1{-# LANGUAGE FlexibleContexts #-}
2{-# LANGUAGE CPP #-}
3-----------------------------------------------------------------------------
4-- |
5-- Module : Data.Packed.Vector
6-- Copyright : (c) Alberto Ruiz 2007-10
7-- License : BSD3
8-- Maintainer : Alberto Ruiz
9-- Stability : provisional
10--
11-- 1D arrays suitable for numeric computations using external libraries.
12--
13-- This module provides basic functions for manipulation of structure.
14--
15-----------------------------------------------------------------------------
16{-# OPTIONS_HADDOCK hide #-}
17
18module Data.Packed.Vector (
19 Vector,
20 fromList, (|>), toList, buildVector,
21 dim, (@>),
22 subVector, takesV, vjoin, join,
23 mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith,
24 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
25 foldLoop, foldVector, foldVectorG, foldVectorWithIndex,
26 toByteString, fromByteString
27) where
28
29import Data.Packed.Internal.Vector
30import Foreign.Storable
31
32-------------------------------------------------------------------
33
34#ifdef BINARY
35
36import Data.Binary
37import Control.Monad(replicateM)
38
39import Data.ByteString.Internal as BS
40import Foreign.ForeignPtr(castForeignPtr)
41import Data.Vector.Storable.Internal(updPtr)
42import Foreign.Ptr(plusPtr)
43
44
45-- a 64K cache, with a Double taking 13 bytes in Bytestring,
46-- implies a chunk size of 5041
47chunk :: Int
48chunk = 5000
49
50chunks :: Int -> [Int]
51chunks d = let c = d `div` chunk
52 m = d `mod` chunk
53 in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk)
54
55putVector v = mapM_ put $! toList v
56
57getVector d = do
58 xs <- replicateM d get
59 return $! fromList xs
60
61--------------------------------------------------------------------------------
62
63toByteString :: Storable t => Vector t -> ByteString
64toByteString v = BS.PS (castForeignPtr fp) (sz*o) (sz * dim v)
65 where
66 (fp,o,_n) = unsafeToForeignPtr v
67 sz = sizeOf (v@>0)
68
69
70fromByteString :: Storable t => ByteString -> Vector t
71fromByteString (BS.PS fp o n) = r
72 where
73 r = unsafeFromForeignPtr (castForeignPtr (updPtr (`plusPtr` o) fp)) 0 n'
74 n' = n `div` sz
75 sz = sizeOf (r@>0)
76
77--------------------------------------------------------------------------------
78
79instance (Binary a, Storable a) => Binary (Vector a) where
80
81 put v = do
82 let d = dim v
83 put d
84 mapM_ putVector $! takesV (chunks d) v
85
86 -- put = put . v2bs
87
88 get = do
89 d <- get
90 vs <- mapM getVector $ chunks d
91 return $! vjoin vs
92
93 -- get = fmap bs2v get
94
95#endif
96
97
98-------------------------------------------------------------------
99
100{- | creates a Vector of the specified length using the supplied function to
101 to map the index to the value at that index.
102
103@> buildVector 4 fromIntegral
1044 |> [0.0,1.0,2.0,3.0]@
105
106-}
107buildVector :: Storable a => Int -> (Int -> a) -> Vector a
108buildVector len f =
109 fromList $ map f [0 .. (len - 1)]
110
111
112-- | zip for Vectors
113zipVector :: (Storable a, Storable b, Storable (a,b)) => Vector a -> Vector b -> Vector (a,b)
114zipVector = zipVectorWith (,)
115
116-- | unzip for Vectors
117unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vector a,Vector b)
118unzipVector = unzipVectorWith id
119
120-------------------------------------------------------------------
121
122{-# DEPRECATED join "use vjoin or Data.Vector.concat" #-}
123join :: Storable t => [Vector t] -> Vector t
124join = vjoin
125
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-20 01:26:06 +0000