Merge remote-tracking branch 'origin/master'

5 files changed, 149 insertions, 95 deletions

diff --git a/lib/Data/Packed/Foreign.hs b/lib/Data/Packed/Foreign.hs
new file mode 100644
index 0000000..efa51ca
--- /dev/null
+++ b/lib/Data/Packed/Foreign.hs
@@ -0,0 +1,99 @@
+{-# LANGUAGE MagicHash, UnboxedTuples #-}
+-- | FFI and hmatrix helpers.
+--
+-- Sample usage, to upload a perspective matrix to a shader.
+--
+-- @ glUniformMatrix4fv 0 1 (fromIntegral gl_TRUE) \`appMatrix\` perspective 0.01 100 (pi\/2) (4\/3) 
+-- @
+--
+module Data.Packed.Foreign 
+    ( app
+    , appVector, appVectorLen
+    , appMatrix, appMatrixLen, appMatrixRaw, appMatrixRawLen
+    , unsafeMatrixToVector, unsafeMatrixToForeignPtr
+    ) where
+import Data.Packed.Internal
+import qualified Data.Vector.Storable as S
+import Foreign (Ptr, ForeignPtr, Storable)
+import Foreign.C.Types (CInt)
+import GHC.Base (IO(..), realWorld#)
+
+{-# INLINE unsafeInlinePerformIO #-}
+-- | 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
+-- unecessary calls to unsafePerformIO or its internals.
+unsafeInlinePerformIO :: IO a -> a
+unsafeInlinePerformIO (IO f) = case f realWorld# of
+    (# _, x #) -> x
+
+{-# INLINE app #-}
+-- | Only useful since it is left associated with a precedence of 1, unlike 'Prelude.$', which is right associative.
+-- e.g.
+--
+-- @
+-- someFunction
+--     \`appMatrixLen\` m
+--     \`appVectorLen\` v
+--     \`app\` other
+--     \`app\` arguments
+--     \`app\` go here
+-- @
+--
+-- One could also write:
+--
+-- @
+-- (someFunction 
+--     \`appMatrixLen\` m
+--     \`appVectorLen\` v) 
+--     other 
+--     arguments 
+--     (go here)
+-- @
+--
+app :: (a -> b) -> a -> b
+app f = f
+
+{-# INLINE appVector #-}
+appVector :: Storable a => (Ptr a -> b) -> Vector a -> b
+appVector f x = unsafeInlinePerformIO (S.unsafeWith x (return . f))
+
+{-# INLINE appVectorLen #-}
+appVectorLen :: Storable a => (CInt -> Ptr a -> b) -> Vector a -> b
+appVectorLen f x = unsafeInlinePerformIO (S.unsafeWith x (return . f (fromIntegral (S.length x))))
+
+{-# INLINE appMatrix #-}
+appMatrix :: Element a => (Ptr a -> b) -> Matrix a -> b
+appMatrix f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f))
+
+{-# INLINE appMatrixLen #-}
+appMatrixLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
+appMatrixLen f x = unsafeInlinePerformIO (S.unsafeWith (flatten x) (return . f r c))
+  where
+    r = fromIntegral (rows x)
+    c = fromIntegral (cols x)
+
+{-# INLINE appMatrixRaw #-}
+appMatrixRaw :: Storable a => (Ptr a -> b) -> Matrix a -> b
+appMatrixRaw f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f))
+
+{-# INLINE appMatrixRawLen #-}
+appMatrixRawLen :: Element a => (CInt -> CInt -> Ptr a -> b) -> Matrix a -> b
+appMatrixRawLen f x = unsafeInlinePerformIO (S.unsafeWith (xdat x) (return . f r c))
+  where
+    r = fromIntegral (rows x)
+    c = fromIntegral (cols x)
+
+infixl 1 `app`
+infixl 1 `appVector`
+infixl 1 `appMatrix`
+infixl 1 `appMatrixRaw`
+
+{-# INLINE unsafeMatrixToVector #-}
+-- | This will disregard the order of the matrix, and simply return it as-is. 
+-- If the order of the matrix is RowMajor, this function is identical to 'flatten'.
+unsafeMatrixToVector :: Matrix a -> Vector a
+unsafeMatrixToVector = xdat
+
+{-# INLINE unsafeMatrixToForeignPtr #-}
+unsafeMatrixToForeignPtr :: Storable a => Matrix a -> (ForeignPtr a, Int)
+unsafeMatrixToForeignPtr m = S.unsafeToForeignPtr0 (xdat m)
+
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index b8ed18d..255009c 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -47,6 +47,7 @@ import Data.Complex(Complex)
 import Foreign.C.Types
 import Foreign.C.String(newCString)
 import System.IO.Unsafe(unsafePerformIO)
+import Control.DeepSeq
 
 -----------------------------------------------------------------
 
@@ -131,12 +132,10 @@ mat a f =
         let m g = do
             g (fi (rows a)) (fi (cols a)) p
         f m
-
-{- | Creates a vector by concatenation of rows
-
-@\> flatten ('ident' 3)
-9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
--}
+-- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
+--
+-- @\> flatten ('ident' 3)
+-- 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@
 flatten :: Element t => Matrix t -> Vector t
 flatten = xdat . cmat
 
@@ -459,3 +458,13 @@ size m = (rows m, cols m)
 
 shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
 
+----------------------------------------------------------------------
+
+instance (Storable t, NFData t) => NFData (Matrix t)
+  where
+    rnf m | d > 0     = rnf (v @> 0)
+          | otherwise = ()
+      where
+        d = dim v
+        v = xdat m
+
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index d3b80ff..5892e67 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -55,48 +55,17 @@ import GHC.Base
 import GHC.IOBase hiding (liftIO)
 #endif
 
-#ifdef VECTOR
 import qualified Data.Vector.Storable as Vector
 import Data.Vector.Storable(Vector,
                             unsafeToForeignPtr,
                             unsafeFromForeignPtr,
                             unsafeWith)
-#else
-import Foreign.ForeignPtr(withForeignPtr)
-#endif
 
-#ifdef VECTOR
 
 -- | Number of elements
 dim :: (Storable t) => Vector t -> Int
 dim = Vector.length
 
-#else
-
--- | One-dimensional array of objects stored in a contiguous memory block.
-data Vector t =
-    V { ioff :: {-# UNPACK #-} !Int              -- ^ offset of first element
-      , idim :: {-# UNPACK #-} !Int              -- ^ number of elements
-      , fptr :: {-# UNPACK #-} !(ForeignPtr t)   -- ^ foreign pointer to the memory block
-      }
-
-unsafeToForeignPtr :: Storable a => Vector a -> (ForeignPtr a, Int, Int)
-unsafeToForeignPtr v = (fptr v, ioff v, idim v)
-
--- | Same convention as in Roman Leshchinskiy's vector package.
-unsafeFromForeignPtr :: Storable a => ForeignPtr a -> Int -> Int -> Vector a
-unsafeFromForeignPtr fp i n | n > 0 = V {ioff = i, idim = n, fptr = fp}
-                            | otherwise = error "unsafeFromForeignPtr with dim < 1"
-
-unsafeWith (V i _ fp) m = withForeignPtr fp $ \p -> m (p `advancePtr` i)
-{-# INLINE unsafeWith #-}
-
--- | Number of elements
-dim :: (Storable t) => Vector t -> Int
-dim = idim
-
-#endif
-
 
 -- C-Haskell vector adapter
 -- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
@@ -204,36 +173,8 @@ subVector :: Storable t => Int       -- ^ index of the starting element
                         -> Int       -- ^ number of elements to extract
                         -> Vector t  -- ^ source
                         -> Vector t  -- ^ result
-
-#ifdef VECTOR
-
 subVector = Vector.slice
 
-{-
-subVector k l v
-    | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
-    | otherwise = unsafeFromForeignPtr fp (i+k) l
-  where
-    (fp, i, n) = unsafeToForeignPtr v
--}
-
-#else
-
-subVector k l v@V{idim = n, ioff = i}
-    | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
-    | otherwise = v {idim = l, ioff = i+k}
-
-{-
-subVectorCopy k l (v@V {idim=n})
-    | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
-    | otherwise = unsafePerformIO $ do
-        r <- createVector l
-        let f _ s _ d = copyArray d (advancePtr s k) l >> return 0
-        app2 f vec v vec r "subVector"
-        return r
--}
-
-#endif
 
 {- | Reads a vector position:
 
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index c1a9b24..1b67820 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -30,7 +30,7 @@ module Data.Packed.Matrix (
     (@@>),
     asRow, asColumn,
     fromRows, toRows, fromColumns, toColumns,
-    fromBlocks, toBlocks, toBlocksEvery,
+    fromBlocks, diagBlock, toBlocks, toBlocksEvery,
     repmat,
     flipud, fliprl,
     subMatrix, takeRows, dropRows, takeColumns, dropColumns,
@@ -46,6 +46,7 @@ import Data.Array
 
 import Data.List(transpose,intersperse)
 import Foreign.Storable(Storable)
+import Control.Monad(liftM)
 
 -------------------------------------------------------------------
 
@@ -155,7 +156,19 @@ adaptBlocks ms = ms' where
         x = m@@>(0,0)
     g _ _ = error "inconsistent dimensions in fromBlocks"
 
------------------------------------------------------------
+
+--------------------------------------------------------------------------------
+
+-- | create a block diagonal matrix
+diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
+diagBlock ms = fromBlocks $ zipWith f ms [0..]
+  where
+    f m k = take n $ replicate k z ++ m : repeat z
+    n = length ms
+    z = (1><1) [0]
+
+--------------------------------------------------------------------------------
+
 
 -- | Reverse rows
 flipud :: Element t => Matrix t -> Matrix t
@@ -351,7 +364,11 @@ toBlocksEvery r c m = toBlocks rs cs m where
 
 -------------------------------------------------------------------
 
-mk c g = \k v -> g (divMod k c) v
+-- Given a column number and a function taking matrix indexes, returns
+-- a function which takes vector indexes (that can be used on the
+-- flattened matrix).
+mk :: Int -> ((Int, Int) -> t) -> (Int -> t)
+mk c g = \k -> g (divMod k c)
 
 {- |
 
@@ -364,9 +381,8 @@ m[1,1] = 5
 m[1,2] = 6@
 -}
 mapMatrixWithIndexM_
-  :: (Element a, Num a,
-      Functor f, Monad f) =>
-      ((Int, Int) -> a -> f ()) -> Matrix a -> f ()
+  :: (Element a, Num a, Monad m) =>
+      ((Int, Int) -> a -> m ()) -> Matrix a -> m ()
 mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
   where
     c = cols m
@@ -380,11 +396,9 @@ Just (3><3)
  ,  20.0,  21.0, 122.0 ]@
 -}
 mapMatrixWithIndexM
-  :: (Foreign.Storable.Storable t,
-      Element a, Num a,
-      Functor f, Monad f) =>
-      ((Int, Int) -> a -> f t) -> Matrix a -> f (Matrix t)
-mapMatrixWithIndexM g m = fmap (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m
+  :: (Element a, Storable b, Monad m) =>
+      ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
+mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m 
     where
       c = cols m
 
@@ -395,10 +409,10 @@ mapMatrixWithIndexM g m = fmap (reshape c) . mapVectorWithIndexM (mk c g) . flat
  ,  10.0, 111.0,  12.0
  ,  20.0,  21.0, 122.0 ]@
  -}
-mapMatrixWithIndex :: (Foreign.Storable.Storable t,
-      Element a, Num a) =>
-      ((Int, Int) -> a -> t) -> Matrix a -> Matrix t
-mapMatrixWithIndex g m = reshape c $ mapVectorWithIndex (mk c g) $ flatten m
+mapMatrixWithIndex
+  :: (Element a, Storable b) =>
+      ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
+mapMatrixWithIndex g m = reshape c . mapVectorWithIndex (mk c g) . flatten $ m
     where
       c = cols m
 
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
index 4b229f0..dabb17d 100644
--- a/lib/Data/Packed/Random.hs
+++ b/lib/Data/Packed/Random.hs
@@ -12,11 +12,11 @@
 -----------------------------------------------------------------------------
 
 module Data.Packed.Random (
+    Seed,
     RandDist(..),
     randomVector,
     gaussianSample,
-    uniformSample,
-    meanCov,
+    uniformSample
 ) where
 
 import Numeric.GSL.Vector
@@ -25,9 +25,11 @@ import Numeric.ContainerBoot
 import Numeric.LinearAlgebra.Algorithms
 
 
+type Seed = Int
+
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- Gaussian distribution.
-gaussianSample :: Int -- ^ seed
+gaussianSample :: Seed
                -> Int -- ^ number of rows
                -> Vector Double -- ^ mean vector
                -> Matrix Double -- ^ covariance matrix
@@ -40,7 +42,7 @@ gaussianSample seed n med cov = m where
 
 -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
 -- uniform distribution.
-uniformSample :: Int -- ^ seed
+uniformSample :: Seed
                -> Int -- ^ number of rows
                -> [(Double,Double)] -- ^ ranges for each column
                -> Matrix Double -- ^ result
@@ -53,14 +55,3 @@ uniformSample seed n rgs = m where
     am = konst 1 n `outer` a
     m = fromColumns (zipWith scale cs dat) `add` am
 
------------- utilities -------------------------------
-
--- | Compute mean vector and covariance matrix of the rows of a matrix.
-meanCov :: Matrix Double -> (Vector Double, Matrix Double)
-meanCov x = (med,cov) where
-    r    = rows x
-    k    = 1 / fromIntegral r
-    med  = konst k r `vXm` x
-    meds = konst 1 r `outer` med
-    xc   = x `sub` meds
-    cov  = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)