move util

1 files changed, 568 insertions, 0 deletions

diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
new file mode 100644
index 0000000..9900770
--- /dev/null
+++ b/packages/base/src/Internal/Util.hs
@@ -0,0 +1,568 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE ViewPatterns #-}
+
+
+-----------------------------------------------------------------------------
+{- |
+Module      :  Internal.Util
+Copyright   :  (c) Alberto Ruiz 2013
+License     :  BSD3
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+-}
+-----------------------------------------------------------------------------
+
+module Internal.Util(
+
+    -- * Convenience functions
+    vector, matrix,
+    disp,
+    formatSparse,
+    approxInt,
+    dispDots,
+    dispBlanks,
+    formatShort,
+    dispShort,
+    zeros, ones,
+    diagl,
+    row,
+    col,
+    (&), (¦), (|||), (——), (===), (#),
+    (?), (¿),
+    Indexable(..), size,
+    Numeric,
+    rand, randn,
+    cross,
+    norm,
+    ℕ,ℤ,ℝ,ℂ,iC,
+    Normed(..), norm_Frob, norm_nuclear,
+    unitary,
+    mt,
+    (~!~),
+    pairwiseD2,
+    rowOuters,
+    null1,
+    null1sym,
+    -- * Convolution
+    -- ** 1D
+    corr, conv, corrMin,
+    -- ** 2D
+    corr2, conv2, separable,
+    gaussElim
+) where
+
+import Internal.Tools
+import Internal.Vector
+import Internal.Matrix hiding (size)
+import Internal.Numeric
+import Internal.Element
+import Internal.Container
+import Internal.Vectorized
+import Internal.IO
+import Internal.Algorithms hiding (i,Normed)
+import Numeric.Matrix()
+import Numeric.Vector()
+import Internal.Random
+import Internal.Convolution
+import Control.Monad(when)
+import Text.Printf
+import Data.List.Split(splitOn)
+import Data.List(intercalate,sortBy)
+import Data.Function(on)
+import Control.Arrow((&&&))
+import Data.Complex
+import Data.Vector.Storable(fromList)
+
+type ℝ = Double
+type ℕ = Int
+type ℤ = Int
+type ℂ = Complex Double
+
+-- | imaginary unit
+iC :: ℂ
+iC = 0:+1
+
+{- | Create a real vector.
+
+>>> vector [1..5]
+fromList [1.0,2.0,3.0,4.0,5.0]
+
+-}
+vector :: [ℝ] -> Vector ℝ
+vector = fromList
+
+{- | Create a real matrix.
+
+>>> matrix 5 [1..15]
+(3><5)
+ [  1.0,  2.0,  3.0,  4.0,  5.0
+ ,  6.0,  7.0,  8.0,  9.0, 10.0
+ , 11.0, 12.0, 13.0, 14.0, 15.0 ]
+
+-}
+matrix
+  :: Int -- ^ number of columns
+  -> [ℝ] -- ^ elements in row order
+  -> Matrix ℝ
+matrix c = reshape c . fromList
+
+
+{- | print a real matrix with given number of digits after the decimal point
+
+>>> disp 5 $ ident 2 / 3
+2x2
+0.33333  0.00000
+0.00000  0.33333
+
+-}
+disp :: Int -> Matrix Double -> IO ()
+
+disp n = putStr . dispf n
+
+
+{- | create a real diagonal matrix from a list
+
+>>> diagl [1,2,3]
+(3><3)
+ [ 1.0, 0.0, 0.0
+ , 0.0, 2.0, 0.0
+ , 0.0, 0.0, 3.0 ]
+
+-}
+diagl :: [Double] -> Matrix Double
+diagl = diag . fromList
+
+-- | a real matrix of zeros
+zeros :: Int -- ^ rows
+      -> Int -- ^ columns
+      -> Matrix Double
+zeros r c = konst 0 (r,c)
+
+-- | a real matrix of ones
+ones :: Int -- ^ rows
+     -> Int -- ^ columns
+     -> Matrix Double
+ones r c = konst 1 (r,c)
+
+-- | concatenation of real vectors
+infixl 3 &
+(&) :: Vector Double -> Vector Double -> Vector Double
+a & b = vjoin [a,b]
+
+{- | horizontal concatenation of real matrices
+
+>>> ident 3 ||| konst 7 (3,4)
+(3><7)
+ [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
+ , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
+ , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]
+
+-}
+infixl 3 |||
+(|||) :: Matrix Double -> Matrix Double -> Matrix Double
+a ||| b = fromBlocks [[a,b]]
+
+-- | a synonym for ('|||') (unicode 0x00a6, broken bar)
+infixl 3 ¦
+(¦) :: Matrix Double -> Matrix Double -> Matrix Double
+(¦) = (|||)
+
+
+-- | vertical concatenation of real matrices
+--
+(===) :: Matrix Double -> Matrix Double -> Matrix Double
+infixl 2 ===
+a === b = fromBlocks [[a],[b]]
+
+-- | a synonym for ('===') (unicode 0x2014, em dash)
+(——) :: Matrix Double -> Matrix Double -> Matrix Double
+infixl 2 ——
+(——) = (===)
+
+
+(#) :: Matrix Double -> Matrix Double -> Matrix Double
+infixl 2 #
+a # b = fromBlocks [[a],[b]]
+
+-- | create a single row real matrix from a list
+--
+-- >>> row [2,3,1,8]
+-- (1><4)
+--  [ 2.0, 3.0, 1.0, 8.0 ]
+--
+row :: [Double] -> Matrix Double
+row = asRow . fromList
+
+-- | create a single column real matrix from a list
+--
+-- >>> col [7,-2,4]
+-- (3><1)
+--  [  7.0
+--  , -2.0
+--  ,  4.0 ]
+--
+col :: [Double] -> Matrix Double
+col = asColumn . fromList
+
+{- | extract rows
+
+>>> (20><4) [1..] ? [2,1,1]
+(3><4)
+ [ 9.0, 10.0, 11.0, 12.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 5.0,  6.0,  7.0,  8.0 ]
+
+-}
+infixl 9 ?
+(?) :: Element t => Matrix t -> [Int] -> Matrix t
+(?) = flip extractRows
+
+{- | extract columns
+
+(unicode 0x00bf, inverted question mark, Alt-Gr ?)
+
+>>> (3><4) [1..] ¿ [3,0]
+(3><2)
+ [  4.0, 1.0
+ ,  8.0, 5.0
+ , 12.0, 9.0 ]
+
+-}
+infixl 9 ¿
+(¿) :: Element t => Matrix t -> [Int] -> Matrix t
+(¿)= flip extractColumns
+
+
+cross :: Product t => Vector t -> Vector t -> Vector t
+-- ^ cross product (for three-element vectors)
+cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
+          | otherwise = error $ "the cross product requires 3-element vectors (sizes given: "
+                                ++show (dim x)++" and "++show (dim y)++")"
+  where
+    [x1,x2,x3] = toList x
+    [y1,y2,y3] = toList y
+    z1 = x2*y3-x3*y2
+    z2 = x3*y1-x1*y3
+    z3 = x1*y2-x2*y1
+
+{-# SPECIALIZE cross :: Vector Double -> Vector Double -> Vector Double #-}
+{-# SPECIALIZE cross :: Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) #-}
+
+norm :: Vector Double -> Double
+-- ^ 2-norm of real vector
+norm = pnorm PNorm2
+
+class Normed a
+  where
+    norm_0   :: a -> ℝ
+    norm_1   :: a -> ℝ
+    norm_2   :: a -> ℝ
+    norm_Inf :: a -> ℝ
+
+
+instance Normed (Vector ℝ)
+  where
+    norm_0 v = sumElements (step (abs v - scalar (eps*normInf v)))
+    norm_1 = pnorm PNorm1
+    norm_2 = pnorm PNorm2
+    norm_Inf = pnorm Infinity
+
+instance Normed (Vector ℂ)
+  where
+    norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v)))
+    norm_1 = pnorm PNorm1
+    norm_2 = pnorm PNorm2
+    norm_Inf = pnorm Infinity
+
+instance Normed (Matrix ℝ)
+  where
+    norm_0 = norm_0 . flatten
+    norm_1 = pnorm PNorm1
+    norm_2 = pnorm PNorm2
+    norm_Inf = pnorm Infinity
+
+instance Normed (Matrix ℂ)
+  where
+    norm_0 = norm_0 . flatten
+    norm_1 = pnorm PNorm1
+    norm_2 = pnorm PNorm2
+    norm_Inf = pnorm Infinity
+
+instance Normed (Vector I)
+  where
+    norm_0 = fromIntegral . sumElements . step . abs
+    norm_1 = fromIntegral . norm1
+    norm_2 v = sqrt . fromIntegral $ dot v v
+    norm_Inf = fromIntegral . normInf
+
+
+
+norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ
+norm_Frob = norm_2 . flatten
+
+norm_nuclear :: Field t => Matrix t -> ℝ
+norm_nuclear = sumElements . singularValues
+
+
+-- | Obtains a vector in the same direction with 2-norm=1
+unitary :: Vector Double -> Vector Double
+unitary v = v / scalar (norm v)
+
+
+-- | trans . inv
+mt :: Matrix Double -> Matrix Double
+mt = trans . inv
+
+--------------------------------------------------------------------------------
+{- |
+
+>>> size $ vector [1..10]
+10
+>>> size $ (2><5)[1..10::Double]
+(2,5)
+
+-}
+size :: Container c t => c t -> IndexOf c
+size = size'
+
+{- | Alternative indexing function.
+
+>>> vector [1..10] ! 3
+4.0
+
+On a matrix it gets the k-th row as a vector:
+
+>>> matrix 5 [1..15] ! 1
+fromList [6.0,7.0,8.0,9.0,10.0]
+
+>>> matrix 5 [1..15] ! 1 ! 3
+9.0
+
+-}
+class Indexable c t | c -> t , t -> c
+  where
+    infixl 9 !
+    (!) :: c -> Int -> t
+
+instance Indexable (Vector Double) Double
+  where
+    (!) = (@>)
+
+instance Indexable (Vector Float) Float
+  where
+    (!) = (@>)
+
+instance Indexable (Vector I) I
+  where
+    (!) = (@>)
+
+instance Indexable (Vector (Complex Double)) (Complex Double)
+  where
+    (!) = (@>)
+
+instance Indexable (Vector (Complex Float)) (Complex Float)
+  where
+    (!) = (@>)
+
+instance Element t => Indexable (Matrix t) (Vector t)
+  where
+    m!j = subVector (j*c) c (flatten m)
+      where
+        c = cols m
+
+--------------------------------------------------------------------------------
+
+-- | Matrix of pairwise squared distances of row vectors
+-- (using the matrix product trick in blog.smola.org)
+pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
+pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
+               | otherwise = error $ "pairwiseD2 with different number of columns: "
+                                   ++ show (size x) ++ ", " ++ show (size y)
+  where
+    ox = one (rows x)
+    oy = one (rows y)
+    oc = one (cols x)
+    one k = konst 1 k
+    x2 = x * x <> oc
+    y2 = y * y <> oc
+    ok = cols x == cols y
+
+--------------------------------------------------------------------------------
+
+{- | outer products of rows
+
+>>> a
+(3><2)
+ [   1.0,   2.0
+ ,  10.0,  20.0
+ , 100.0, 200.0 ]
+>>> b
+(3><3)
+ [ 1.0, 2.0, 3.0
+ , 4.0, 5.0, 6.0
+ , 7.0, 8.0, 9.0 ]
+
+>>> rowOuters a (b ||| 1)
+(3><8)
+ [   1.0,   2.0,   3.0,   1.0,    2.0,    4.0,    6.0,   2.0
+ ,  40.0,  50.0,  60.0,  10.0,   80.0,  100.0,  120.0,  20.0
+ , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ]
+
+-}
+rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
+rowOuters a b = a' * b'
+  where
+    a' = kronecker a (ones 1 (cols b))
+    b' = kronecker (ones 1 (cols a)) b
+
+--------------------------------------------------------------------------------
+
+-- | solution of overconstrained homogeneous linear system
+null1 :: Matrix Double -> Vector Double
+null1 = last . toColumns . snd . rightSV
+
+-- | solution of overconstrained homogeneous symmetric linear system
+null1sym :: Matrix Double -> Vector Double
+null1sym = last . toColumns . snd . eigSH'
+
+--------------------------------------------------------------------------------
+
+infixl 0 ~!~
+c ~!~ msg = when c (error msg)
+
+--------------------------------------------------------------------------------
+
+formatSparse :: String -> String -> String -> Int -> Matrix Double -> String
+
+formatSparse zeroI _zeroF sep _ (approxInt -> Just m) = format sep f m
+  where
+    f 0 = zeroI
+    f x = printf "%.0f" x
+
+formatSparse zeroI zeroF sep n m = format sep f m
+  where
+    f x | abs (x::Double) < 2*peps = zeroI++zeroF
+        | abs (fromIntegral (round x::Int) - x) / abs x < 2*peps
+            = printf ("%.0f."++replicate n ' ') x
+        | otherwise = printf ("%."++show n++"f") x
+
+approxInt m
+    | norm_Inf (v - vi) < 2*peps * norm_Inf v = Just (reshape (cols m) vi)
+    | otherwise = Nothing
+  where
+    v = flatten m
+    vi = roundVector v
+
+dispDots n = putStr . formatSparse "." (replicate n ' ') "  " n
+
+dispBlanks n = putStr . formatSparse "" "" "  " n
+
+formatShort sep fmt maxr maxc m = auxm4
+  where
+    (rm,cm) = size m
+    (r1,r2,r3)
+        | rm <= maxr = (rm,0,0)
+        | otherwise  = (maxr-3,rm-maxr+1,2)
+    (c1,c2,c3)
+        | cm <= maxc = (cm,0,0)
+        | otherwise  = (maxc-3,cm-maxc+1,2)
+    [ [a,_,b]
+     ,[_,_,_]
+     ,[c,_,d]] = toBlocks [r1,r2,r3]
+                          [c1,c2,c3] m
+    auxm = fromBlocks [[a,b],[c,d]]
+    auxm2
+        | cm > maxc = format "|" fmt auxm
+        | otherwise = format sep fmt auxm
+    auxm3
+        | cm > maxc = map (f . splitOn "|") (lines auxm2)
+        | otherwise = (lines auxm2)
+    f items = intercalate sep (take (maxc-3) items) ++ "  .. " ++
+              intercalate sep (drop (maxc-3) items)
+    auxm4
+        | rm > maxr = unlines (take (maxr-3) auxm3 ++ vsep : drop (maxr-3) auxm3)
+        | otherwise = unlines auxm3
+    vsep = map g (head auxm3)
+    g '.' = ':'
+    g _ = ' '
+
+
+dispShort :: Int -> Int -> Int -> Matrix Double -> IO ()
+dispShort maxr maxc dec m =
+    printf "%dx%d\n%s" (rows m) (cols m) (formatShort "  " fmt maxr maxc m)
+  where
+    fmt = printf ("%."++show dec ++"f")
+
+--------------------------------------------------------------------------------
+
+-- | generic reference implementation of gaussian elimination
+--
+-- @a <> gauss a b = b@
+--
+gaussElim
+  :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
+  => Matrix t -> Matrix t -> Matrix t
+
+gaussElim x y = dropColumns (rows x) (flipud $ fromRows s2)
+  where
+    rs = toRows $ fromBlocks [[x , y]]
+    s1 = pivotDown (rows x) 0 rs
+    s2 = pivotUp (rows x-1) (reverse s1)
+
+pivotDown t n xs
+    | t == n    = []
+    | otherwise = y : pivotDown t (n+1) ys
+  where
+    y:ys = redu (pivot n xs)
+
+    pivot k = (const k &&& id)
+            . reverse . sortBy (compare `on` (abs. (!k)))  -- FIXME
+
+    redu (k,x:zs)
+        | p == 0 = error "gauss: singular!"  -- FIXME
+        | otherwise = u : map f zs
+      where
+        p = x!k
+        u = scale (recip (x!k)) x
+        f z = z - scale (z!k) u
+    redu (_,[]) = []
+
+
+pivotUp n xs
+    | n == -1 = []
+    | otherwise = y : pivotUp (n-1) ys
+  where
+    y:ys = redu' (n,xs)
+
+    redu' (k,x:zs) = u : map f zs
+      where
+        u = x
+        f z = z - scale (z!k) u
+    redu' (_,[]) = []
+
+--------------------------------------------------------------------------------
+
+instance Testable (Matrix I) where
+   checkT _ = test
+
+test :: (Bool, IO())
+test = (and ok, return ())
+  where
+    m  = (3><4) [1..12] :: Matrix I
+    r  = (2><3) [1,2,3,4,3,2]
+    c  = (3><2) [0,4,4,1,2,3]
+    p  = (9><10) [0..89] :: Matrix I
+    ep = (2><3) [10,24,32,44,31,23]
+    md = fromInt m      :: Matrix Double
+    ok = [ tr m <> m == toInt (tr md <> md)
+         , m <> tr m == toInt (md <> tr md)
+         , m ?? (Take 2, Take 3) == remap (asColumn (range 2)) (asRow (range 3)) m
+         , remap r (tr c) p == ep
+         , tr p ?? (PosCyc (idxs[-5,13]), Pos (idxs[3,7,1])) == (2><3) [35,75,15,33,73,13]
+         ]
+