7 files changed, 22 insertions, 312 deletions

diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
deleted file mode 100644
index d7f8c40..0000000
--- a/examples/benchmarks.hs
+++ /dev/null
@@ -1,151 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-
--- $ ghc --make -O2 benchmarks.hs
-
-
-import Numeric.LinearAlgebra
-import System.Time
-import System.CPUTime
-import Text.Printf
-import Data.List(foldl1')
-
-
-time act = do
-    t0 <- getCPUTime
-    act
-    t1 <- getCPUTime
-    printf "%.3f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()
-
---------------------------------------------------------------------------------
-
-main = sequence_ [bench1,
-                  bench2,
-                  bench4,
-                  bench5 1000000 3, bench5 100000 50,
-                  bench6 100 (100000::Double), bench6 100000 (100::Double), bench6 10000 (1000::Double)]
-
-w :: Vector Double
-w = constant 1 5000000
-w2 = 1 * w
-
-v = flatten $ ident 500 :: Vector Double
-
-
-bench1 = do
-    time $ print$ vectorMax (w+w2) -- evaluate it
-    putStrLn "Sum of a vector with 5M doubles:"
-    print$ vectorMax (w+w2) -- evaluate it
-    time $ printf "         BLAS: %.2f: " $ sumVB w
-    time $ printf "BLAS only dot: %.2f: " $ w <.> w2
-    time $ printf "      Haskell: %.2f: " $ sumVH w
-    time $ printf "       innerH: %.2f: " $ innerH w w2
-    time $ printf "foldVector: %.2f: " $ sumVector w
-    let getPos k s = if k `mod` 500 < 200 && w@>k > 0 then k:s else s
-    putStrLn "foldLoop for element selection:"
-    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim w)
-    putStrLn "constant 5M:"
-    time $ print $ constant (1::Double) 5000001 @> 7
-    time $ print $ constant           i 5000001 @> 7
-    time $ print $ conj (constant i 5000001) @> 7
-    putStrLn "zips C vs H:"
-    time $ print $ (w / w2) @> 7
-    time $ print $ (zipVector (/) w w2) @> 7
-    putStrLn "folds C/BLAS vs H:"
-    let t = constant (1::Double) 5000002
-    print $ t @> 7
-    time $ print $ foldVector max (t@>0) t
-    time $ print $ vectorMax t
-    time $ print $ sqrt $ foldVector (\v s -> v*v+s) 0 t
-    time $ print $ pnorm PNorm2 t
-    putStrLn "scale C/BLAS vs H:"
-    time $ print $ mapVector (*2) t @> 7
-    time $ print $ (2 * t) @> 7
-
-sumVB v = constant 1 (dim v) <.> v
-
-sumVH v = go (d - 1) 0
-     where
-       d = dim v
-       go :: Int -> Double -> Double
-       go 0 s = s + (v @> 0)
-       go !j !s = go (j - 1) (s + (v @> j))
-
-innerH u v = go (d - 1) 0
-     where
-       d = min (dim u) (dim v)
-       go :: Int -> Double -> Double
-       go 0 s = s + (u @> 0) * (v @> 0)
-       go !j !s = go (j - 1) (s + (u @> j) * (v @> j))
-
-
--- sumVector = foldVectorG (\k v s -> v k + s) 0.0
-sumVector = foldVector (+) 0.0
-
---------------------------------------------------------------------------------
-
-bench2 = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "Multiplication of 1M different 3x3 matrices:"
---    putStrLn "from [[]]"
---    time $ print $ manymult (10^6) rot'
---    putStrLn "from (3><3) []"
-    time $ print $ manymult (10^6) rot
-    print $ cos (10^6/2)
-
-
-rot' :: Double -> Matrix Double
-rot' a = matrix [[ c,0,s],
-                 [ 0,1,0],
-                 [-s,0,c]]
-    where c = cos a
-          s = sin a
-          matrix = fromLists
-
-rot :: Double -> Matrix Double
-rot a = (3><3) [ c,0,s
-               , 0,1,0
-               ,-s,0,c ]
-    where c = cos a
-          s = sin a
-
-manymult n r = foldl1' (<>) (map r angles)
-    where angles = toList $ linspace n (0,1)
-          -- angles = map (k*) [0..n']
-          -- n' = fromIntegral n - 1
-          -- k  = recip n'
-
---------------------------------------------------------------------------------
-
-bench4 = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "1000x1000 inverse"
-    let a = ident 1000 :: Matrix Double
-    let b = 2*a
-    print $ vectorMax $ flatten (a+b) -- evaluate it
-    time $ print $ vectorMax $ flatten $ linearSolve a b
-
---------------------------------------------------------------------------------
-
-op1 a b = a <> trans b
-op2 a b = a  + trans b
-
-timep = time . print . vectorMax . flatten
-
-bench5 n d = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "transpose in add"
-    let ms = replicate n ((ident d :: Matrix Double))
-    timep $ foldl1' (+)  ms
-    timep $ foldl1' op2  ms
-    putStrLn "-------------------------------------------------------"
-    putStrLn "transpose in multiply"
-
-    timep $ foldl1' (<>) ms
-    timep $ foldl1' op1  ms
-
---------------------------------------------------------------------------------
-
-bench6 sz n = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "many constants"
-    time $ print $ sum $ map ((@>0). flip constant sz) [1..n]
\ No newline at end of file
diff --git a/examples/latexmat.hs b/examples/latexmat.hs
deleted file mode 100644
index d912a28..0000000
--- a/examples/latexmat.hs
+++ /dev/null
@@ -1,11 +0,0 @@
-import Numeric.LinearAlgebra
-import Text.Printf
-
-disp w l fmt m = unlines $ map (++l) $ lines $ format w (printf fmt) m
-
-latex fmt m = "\\begin{bmatrix}\n" ++ disp " & " " \\\\" fmt m ++ "\\end{bmatrix}"
-
-main = do
-    let m = (3><4) [1..12::Double]
-    putStrLn $ disp " | " "" "%.2f" m
-    putStrLn $ latex "%.3f" m
diff --git a/examples/pca1.hs b/examples/pca1.hs
index 775bb4c..58b5577 100644
--- a/examples/pca1.hs
+++ b/examples/pca1.hs
@@ -8,15 +8,14 @@ import Control.Monad(when)
 type Vec = Vector Double
 type Mat = Matrix Double
 
-sumColumns m = constant 1 (rows m) <> m
 
--- Vec with the mean value of the columns of a Mat
-mean x = sumColumns x / fromIntegral (rows x)
+-- Vector with the mean value of the columns of a matrix
+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
+
+-- covariance matrix of a list of observations stored as rows
+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
+    where xc = x - asRow (mean x)
 
--- covariance Mat of a list of observations as rows of a Mat
-cov x = (trans xc <> xc) / fromIntegral (rows x -1)
-    where xc = center x
-          center m = m - constant 1 (rows m) `outer` mean m
 
 -- creates the compression and decompression functions from the desired number of components
 pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec)
@@ -38,7 +37,7 @@ main = do
         system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
         system("gunzip mnist.txt.gz")
         return ()
-    m <- fromFile "mnist.txt" (5000,785)
+    m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)
     let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)
     let x = toRows xs !! 4  -- an arbitrary test Vec
     let (pe,pd) = pca 10 xs
diff --git a/examples/pca2.hs b/examples/pca2.hs
index 8c20370..c38857c 100644
--- a/examples/pca2.hs
+++ b/examples/pca2.hs
@@ -9,33 +9,31 @@ import Control.Monad(when)
 type Vec = Vector Double
 type Mat = Matrix Double
 
-sumColumns m = constant 1 (rows m) <> m
+-- Vector with the mean value of the columns of a matrix
+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
 
--- Vector with the mean value of the columns of a Mat
-mean x = sumColumns x / fromIntegral (rows x)
+-- covariance matrix of a list of observations stored as rows
+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
+    where xc = x - asRow (mean x)
 
--- covariance matrix of a list of observations as rows of a matrix
-cov x = (trans xc <> xc) / fromIntegral (rows x -1) 
-    where xc = center x
-          center m = m - constant 1 (rows m) `outer` mean m
 
 type Stat = (Vec, [Double], Mat)
--- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)  
+-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)
 stat :: Mat -> Stat
-stat x = (m, toList s, trans v) where   
+stat x = (m, toList s, trans v) where
     m = mean x
-    (s,v) = eigSH' (cov x)   
+    (s,v) = eigSH' (cov x)
 
--- creates the compression and decompression functions from the desired reconstruction 
+-- creates the compression and decompression functions from the desired reconstruction
 -- quality and the statistics of a data set
 pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec)
 pca prec (m,s,v) = (encode,decode)    
   where    
     encode x = vp <> (x - m)
     decode x = x <> vp + m
-    vp = takeRows n v    
+    vp = takeRows n v
     n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)
-    cumSum = tail . scanl (+) 0.0     
+    cumSum = tail . scanl (+) 0.0
     prec' = if prec <=0.0 || prec >= 1.0
                 then error "the precision in pca must be 0<prec<1"
                 else prec
@@ -49,7 +47,7 @@ test st prec x = do
     let (pe,pd) = pca prec st
     let y = pe x
     print $ dim y
-    shdigit (pd y)    
+    shdigit (pd y)
 
 main = do
     ok <- doesFileExist ("mnist.txt")
@@ -58,7 +56,7 @@ main = do
         system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
         system("gunzip mnist.txt.gz")
         return ()
-    m <- fromFile "mnist.txt" (5000,785)
+    m <- loadMatrix "mnist.txt"
     let xs = takeColumns (cols m -1) m
     let x = toRows xs !! 4  -- an arbitrary test vector
     shdigit x
diff --git a/examples/pinv.hs b/examples/pinv.hs
index cab8fc6..7de50b8 100644
--- a/examples/pinv.hs
+++ b/examples/pinv.hs
@@ -3,7 +3,7 @@ import Graphics.Plot
 import Text.Printf(printf)
 
 expand :: Int -> Vector Double -> Matrix Double
-expand n x = fromColumns $ constant 1 (dim x): map (x^) [1 .. n]
+expand n x = fromColumns $ map (x^) [0 .. n]
 
 polynomialModel :: Vector Double -> Vector Double -> Int
                 -> (Vector Double -> Vector Double)
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 49d7b05..0fb772f 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -33,9 +33,7 @@ extra-source-files: examples/tests.hs
                     examples/kalman.hs
                     examples/parallel.hs
                     examples/plot.hs
-                    examples/latexmat.hs
                     examples/inplace.hs
-                    examples/benchmarks.hs
                     examples/error.hs
                     examples/devel/wrappers.hs
                     examples/devel/functions.c
@@ -128,8 +126,7 @@ library
                         Data.Packed.Convert,
                         Data.Packed.ST,
                         Data.Packed.Development,
-                        Data.Packed.Random,
-                        Numeric.LinearAlgebra.Util
+                        Data.Packed.Random
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Signatures,
diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs
deleted file mode 100644
index 24317e3..0000000
--- a/lib/Numeric/LinearAlgebra/Util.hs
+++ /dev/null
@@ -1,122 +0,0 @@
------------------------------------------------------------------------------
-{- |
-Module      :  Numeric.LinearAlgebra.Util
-Copyright   :  (c) Alberto Ruiz 2010
-License     :  GPL
-
-Maintainer  :  Alberto Ruiz (aruiz at um dot es)
-Stability   :  provisional
-Portability :  portable
-
-Alternative interface and utilities for creation of real arrays, useful to work in interactive mode.
-
--}
------------------------------------------------------------------------------
-
-module Numeric.LinearAlgebra.Util(
-    module Numeric.LinearAlgebra,
-    (<>), (*>), (<*), (<\>), (\>),
-    vector,
-    eye,
-    zeros, ones,
-    diagl,
-    row,
-    col,
-    (#),(&), (//), blocks,
-    rand,
-    splitEvery,
-    table
-) where
-
-import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/))
-import Data.Packed.Internal.Common(table,splitEvery)
-import System.Random(randomIO)
-
-infixl 7 <>
--- | Matrix product ('multiply')
-(<>) :: Field t => Matrix t -> Matrix t -> Matrix t
-(<>) = multiply
-
-infixl 7 *>
--- | matrix x vector
-(*>) :: Field t => Matrix t -> Vector t -> Vector t
-m *> v = flatten $ m <> (asColumn v)
-
-infixl 7 <*
--- | vector x matrix
-(<*) :: Field t => Vector t -> Matrix t -> Vector t
-v <* m = flatten $ (asRow v) <> m
-
-
--- | Least squares solution of a linear system for several right-hand sides, similar to the \\ operator of Matlab\/Octave. (\<\\\>) = 'linearSolveSVD'.
-(<\>) :: (Field a) => Matrix a -> Matrix a -> Matrix a
-infixl 7 <\>
-(<\>) = linearSolveSVD
-
--- | Least squares solution of a linear system for a single right-hand side. See '(\<\\\>)'.
-(\>) :: (Field a) => Matrix a -> Vector a -> Vector a
-infixl 7 \>
-m \> v = flatten (m <\> reshape 1 v)
-
--- | Pseudorandom matrix with uniform elements between 0 and 1.
-rand :: Int -- ^ rows
-     -> Int -- ^ columns
-     -> IO (Matrix Double)
-rand r c = do
-    seed <- randomIO
-    return (reshape c $ randomVector seed Uniform (r*c))
-
--- | Real identity matrix.
-eye :: Int -> Matrix Double
-eye = ident
-
--- | Create a real vector from a list.
-vector :: [Double] -> Vector Double
-vector = fromList
-
--- | Create a real diagonal matrix from a list.
-diagl :: [Double] -> Matrix Double
-diagl = diag . vector
-
--- | Create a matrix or zeros.
-zeros :: Int -- ^ rows
-      -> Int -- ^ columns
-      -> Matrix Double
-zeros r c = reshape c (constant 0 (r*c))
-
--- | Create a matrix or ones.
-ones :: Int -- ^ rows
-     -> Int -- ^ columns
-     -> Matrix Double
-ones r c = reshape c (constant 1 (r*c))
-
-
--- | Concatenation of real vectors.
-infixl 9 #
-(#) :: Vector Double -> Vector Double -> Vector Double
-a # b = join [a,b]
-
-
--- | Horizontal concatenation of real matrices.
-infixl 8 &
-(&) :: Matrix Double -> Matrix Double -> Matrix Double
-a & b = fromBlocks [[a,b]]
-
--- | Vertical concatenation of real matrices.
-infixl 7 //
-(//) :: Matrix Double -> Matrix Double -> Matrix Double
-a // b = fromBlocks [[a],[b]]
-
--- | Real block matrix from a rectangular list of lists.
-blocks :: [[Matrix Double]] -> Matrix Double
-blocks = fromBlocks
-
--- | A real matrix with a single row, create from a list of elements.
-row :: [Double] -> Matrix Double
-row = asRow . vector
-
--- | A real matrix with a single column, created from a list of elements.
-col :: [Double] -> Matrix Double
-col = asColumn . vector
-
-