move util

1 files changed, 0 insertions, 560 deletions

diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
deleted file mode 100644
index 37cdc0f..0000000
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ /dev/null
@@ -1,560 +0,0 @@
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FunctionalDependencies #-}
-{-# LANGUAGE ViewPatterns #-}
-
-
------------------------------------------------------------------------------
-{- |
-Module      :  Numeric.LinearAlgebra.Util
-Copyright   :  (c) Alberto Ruiz 2013
-License     :  BSD3
-Maintainer  :  Alberto Ruiz
-Stability   :  provisional
-
--}
------------------------------------------------------------------------------
-
-module Numeric.LinearAlgebra.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 Data.Packed.Numeric
-import Numeric.LinearAlgebra.Algorithms hiding (i,Normed)
---import qualified Numeric.LinearAlgebra.Algorithms as A
-import Numeric.Matrix()
-import Numeric.Vector()
-import Numeric.LinearAlgebra.Random
-import Numeric.LinearAlgebra.Util.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((&&&))
-
-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]
-         ]
-