backward compatibility

16 files changed, 683 insertions, 464 deletions

diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index a92409b..01e3c26 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -38,18 +38,18 @@ library
                         Data.Packed.ST,
                         Data.Packed.Development,
 
+                        Numeric.LinearAlgebra
                         Numeric.LinearAlgebra.LAPACK
                         Numeric.LinearAlgebra.Algorithms
                         Numeric.Container
                         Numeric.LinearAlgebra.Util
 
-                        Numeric.LinearAlgebra
+                        Numeric.HMatrix
                         Numeric.LinearAlgebra.Devel
                         Numeric.LinearAlgebra.Data
 
-    --                  Numeric.LinearAlgebra.Compat
-                        
-                        
+
+
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Signatures,
@@ -61,6 +61,7 @@ library
                         Numeric.Vector
                         Numeric.Matrix
                         Data.Packed.Internal.Numeric
+                        Data.Packed.Numeric
                         Numeric.LinearAlgebra.Util.Convolution
                         Numeric.LinearAlgebra.Util.CG
                         Numeric.LinearAlgebra.Random
@@ -69,7 +70,7 @@ library
 
     C-sources:          src/C/lapack-aux.c
                         src/C/vector-aux.c
-                        
+
 
     extensions:         ForeignFunctionInterface,
                         CPP
diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs
index 1efedc9..72eb16b 100644
--- a/packages/base/src/Data/Packed/Development.hs
+++ b/packages/base/src/Data/Packed/Development.hs
@@ -25,8 +25,7 @@ module Data.Packed.Development (
     unsafeFromForeignPtr,
     unsafeToForeignPtr,
     check, (//),
-    at', atM', fi, table,
-    conformMs, conformVs, shSize, splitEvery
+    at', atM', fi
 ) where
 
 import Data.Packed.Internal
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs
index 94fb30a..f7afa80 100644
--- a/packages/base/src/Data/Packed/IO.hs
+++ b/packages/base/src/Data/Packed/IO.hs
@@ -18,12 +18,12 @@ module Data.Packed.IO (
 ) where
 
 import Data.Packed
-import Data.Packed.Development
 import Text.Printf(printf)
 import Data.List(intersperse)
 import Data.Complex
 import Numeric.Vectorized(vectorScan,saveMatrix)
 import Control.Applicative((<$>))
+import Data.Packed.Internal
 
 {- | Creates a string from a matrix given a separator and a function to show each entry. Using
 this function the user can easily define any desired display function:
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 9cd18df..3c1c1d0 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -49,6 +49,7 @@ import Data.Complex
 import Control.Applicative((<*>))
 
 import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
+import Data.Packed.Internal
 
 -------------------------------------------------------------------
 
@@ -65,7 +66,9 @@ type instance ArgOf Matrix a = a -> a -> a
 -------------------------------------------------------------------
 
 -- | Basic element-by-element functions for numeric containers
-class (Complexable c, Fractional e, Element e) => Container c e where
+class (Complexable c, Fractional e, Element e) => Container c e
+  where
+    size'        :: c e -> IndexOf c
     scalar'      :: e -> c e
     conj'        :: c e -> c e
     scale'       :: e -> c e -> c e
@@ -114,7 +117,9 @@ class (Complexable c, Fractional e, Element e) => Container c e where
 
 --------------------------------------------------------------------------
 
-instance Container Vector Float where
+instance Container Vector Float
+  where
+    size' = dim
     scale' = vectorMapValF Scale
     scaleRecip = vectorMapValF Recip
     addConstant = vectorMapValF AddConstant
@@ -125,7 +130,7 @@ instance Container Vector Float where
     equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
     arctan2' = vectorZipF ATan2
     scalar' x = fromList [x]
-    konst' = constant
+    konst' = constantD
     build' = buildV
     conj' = id
     cmap' = mapVector
@@ -142,7 +147,9 @@ instance Container Vector Float where
     accum' = accumV
     cond' = condV condF
 
-instance Container Vector Double where
+instance Container Vector Double
+  where
+    size' = dim
     scale' = vectorMapValR Scale
     scaleRecip = vectorMapValR Recip
     addConstant = vectorMapValR AddConstant
@@ -153,7 +160,7 @@ instance Container Vector Double where
     equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
     arctan2' = vectorZipR ATan2
     scalar' x = fromList [x]
-    konst' = constant
+    konst' = constantD
     build' = buildV
     conj' = id
     cmap' = mapVector
@@ -170,7 +177,9 @@ instance Container Vector Double where
     accum' = accumV
     cond' = condV condD
 
-instance Container Vector (Complex Double) where
+instance Container Vector (Complex Double)
+  where
+    size' = dim
     scale' = vectorMapValC Scale
     scaleRecip = vectorMapValC Recip
     addConstant = vectorMapValC AddConstant
@@ -181,7 +190,7 @@ instance Container Vector (Complex Double) where
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
     arctan2' = vectorZipC ATan2
     scalar' x = fromList [x]
-    konst' = constant
+    konst' = constantD
     build' = buildV
     conj' = conjugateC
     cmap' = mapVector
@@ -198,7 +207,9 @@ instance Container Vector (Complex Double) where
     accum' = accumV
     cond' = undefined -- cannot match
 
-instance Container Vector (Complex Float) where
+instance Container Vector (Complex Float)
+  where
+    size' = dim
     scale' = vectorMapValQ Scale
     scaleRecip = vectorMapValQ Recip
     addConstant = vectorMapValQ AddConstant
@@ -209,7 +220,7 @@ instance Container Vector (Complex Float) where
     equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
     arctan2' = vectorZipQ ATan2
     scalar' x = fromList [x]
-    konst' = constant
+    konst' = constantD
     build' = buildV
     conj' = conjugateQ
     cmap' = mapVector
@@ -228,7 +239,9 @@ instance Container Vector (Complex Float) where
 
 ---------------------------------------------------------------
 
-instance (Container Vector a) => Container Matrix a where
+instance (Container Vector a) => Container Matrix a
+  where
+    size' = size
     scale' x = liftMatrix (scale' x)
     scaleRecip x = liftMatrix (scaleRecip x)
     addConstant x = liftMatrix (addConstant x)
@@ -637,7 +650,7 @@ diag v = diagRect 0 v n n where n = dim v
 
 -- | creates the identity matrix of given dimension
 ident :: (Num a, Element a) => Int -> Matrix a
-ident n = diag (constant 1 n)
+ident n = diag (constantD 1 n)
 
 --------------------------------------------------------
 
@@ -681,8 +694,12 @@ condV f a b l e t = f a' b' l' e' t'
 
 class Transposable t
   where
+    -- | (conjugate) transpose
     tr :: t -> t
 
+instance (Container Vector t) => Transposable (Matrix t)
+  where
+    tr = ctrans
 
 class Linear t v
   where
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
index b3be823..2acb31a 100644
--- a/packages/base/src/Data/Packed/Matrix.hs
+++ b/packages/base/src/Data/Packed/Matrix.hs
@@ -226,15 +226,12 @@ takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (co
 
 ------------------------------------------------------------
 
-{- | An easy way to create a matrix:
+{- | create a general matrix
 
->>> (2><3)[2,4,7,-3,11,0]
+>>> (2><3) [2, 4, 7+2*𝑖,   -3, 11, 0]
 (2><3)
- [  2.0,  4.0, 7.0
- , -3.0, 11.0, 0.0 ]
-
-This is the format produced by the instances of Show (Matrix a), which
-can also be used for input.
+ [       2.0 :+ 0.0,  4.0 :+ 0.0, 7.0 :+ 2.0
+ , (-3.0) :+ (-0.0), 11.0 :+ 0.0, 0.0 :+ 0.0 ]
 
 The input list is explicitly truncated, so that it can
 safely be used with lists that are too long (like infinite lists).
@@ -244,6 +241,8 @@ safely be used with lists that are too long (like infinite lists).
  [ 1.0, 2.0, 3.0
  , 4.0, 5.0, 6.0 ]
 
+This is the format produced by the instances of Show (Matrix a), which
+can also be used for input.
 
 -}
 (><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
new file mode 100644
index 0000000..01cf6c5
--- /dev/null
+++ b/packages/base/src/Data/Packed/Numeric.hs
@@ -0,0 +1,288 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Numeric
+-- Copyright   :  (c) Alberto Ruiz 2010-14
+-- License     :  BSD3
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
+--
+-- The 'Container' class is used to define optimized generic functions which work
+-- on 'Vector' and 'Matrix' with real or complex elements.
+--
+-- Some of these functions are also available in the instances of the standard
+-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
+--
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Packed.Numeric (
+    -- * Basic functions
+    module Data.Packed,
+    Konst(..), Build(..),
+    linspace,
+    diag, ident,
+    ctrans,
+    -- * Generic operations
+    Container(..),
+    -- add, mul, sub, divide, equal, scaleRecip, addConstant,
+    scalar, conj, scale, arctan2, cmap,
+    atIndex, minIndex, maxIndex, minElement, maxElement,
+    sumElements, prodElements,
+    step, cond, find, assoc, accum,
+    Transposable(..), Linear(..),
+    -- * Matrix product
+    Product(..), udot, dot, (◇),
+    Mul(..),
+    Contraction(..),(<.>),
+    optimiseMult,
+    mXm,mXv,vXm,LSDiv,(<\>),
+    outer, kronecker,
+    -- * Random numbers
+    RandDist(..),
+    randomVector,
+    gaussianSample,
+    uniformSample,
+    meanCov,
+    -- * Element conversion
+    Convert(..),
+    Complexable(),
+    RealElement(),
+
+    RealOf, ComplexOf, SingleOf, DoubleOf,
+
+    IndexOf,
+    module Data.Complex,
+    -- * IO
+    module Data.Packed.IO,
+    -- * Misc
+    Testable(..)
+) where
+
+import Data.Packed
+import Data.Packed.Internal.Numeric
+import Data.Complex
+import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
+import Data.Monoid(Monoid(mconcat))
+import Data.Packed.IO
+import Numeric.LinearAlgebra.Random
+
+------------------------------------------------------------------
+
+{- | Creates a real vector containing a range of values:
+
+>>> linspace 5 (-3,7::Double)
+fromList [-3.0,-0.5,2.0,4.5,7.0]@
+
+>>> linspace 5 (8,2+i) :: Vector (Complex Double)
+fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
+
+Logarithmic spacing can be defined as follows:
+
+@logspace n (a,b) = 10 ** linspace n (a,b)@
+-}
+linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
+linspace 0 (a,b) = fromList[(a+b)/2]
+linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
+    where s = (b-a)/fromIntegral (n-1)
+
+--------------------------------------------------------
+
+{- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction')
+
+(This operator can also be written using the unicode symbol ◇ (25c7).)
+
+Examples:
+
+>>> let a = (3><4)   [1..]      :: Matrix Double
+>>> let v = fromList [1,0,2,-1] :: Vector Double
+>>> let u = fromList [1,2,3]    :: Vector Double
+
+>>> a
+(3><4)
+ [ 1.0,  2.0,  3.0,  4.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 9.0, 10.0, 11.0, 12.0 ]
+
+matrix × matrix:
+
+>>> disp 2 (a <.> trans a)
+3x3
+ 30   70  110
+ 70  174  278
+110  278  446
+
+matrix × vector:
+
+>>> a <.> v
+fromList [3.0,11.0,19.0]
+
+dot product:
+
+>>> u <.> fromList[3,2,1::Double]
+10
+
+For complex vectors the first argument is conjugated:
+
+>>> fromList [1,i] <.> fromList[2*i+1,3]
+1.0 :+ (-1.0)
+
+>>> fromList [1,i,1-i] <.> complex a
+fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
+-}
+infixl 7 <.>
+(<.>) :: Contraction a b c => a -> b -> c
+(<.>) = contraction
+
+
+class Contraction a b c | a b -> c
+  where
+    -- | Matrix product, matrix - vector product, and dot product
+    contraction :: a -> b -> c
+
+instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
+    u `contraction` v = conj u `udot` v
+
+instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
+    contraction = mXv
+
+instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
+    contraction v m = (conj v) `vXm` m
+
+instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
+    contraction = mXm
+
+
+--------------------------------------------------------------------------------
+
+class Mul a b c | a b -> c where
+ infixl 7 <>
+ -- | Matrix-matrix, matrix-vector, and vector-matrix products.
+ (<>)  :: Product t => a t -> b t -> c t
+
+instance Mul Matrix Matrix Matrix where
+    (<>) = mXm
+
+instance Mul Matrix Vector Vector where
+    (<>) m v = flatten $ m <> asColumn v
+
+instance Mul Vector Matrix Vector where
+    (<>) v m = flatten $ asRow v <> m
+
+--------------------------------------------------------------------------------
+
+-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
+infixl 7 <\>
+(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
+(<\>) = linSolve
+
+class LSDiv c
+  where
+    linSolve :: Field t => Matrix t -> c t -> c t
+
+instance LSDiv Vector
+  where
+    linSolve m v = flatten (linearSolveSVD m (reshape 1 v))
+
+instance LSDiv Matrix
+  where
+    linSolve = linearSolveSVD
+
+--------------------------------------------------------------------------------
+
+class Konst e d c | d -> c, c -> d
+  where
+    -- |
+    -- >>> konst 7 3 :: Vector Float
+    -- fromList [7.0,7.0,7.0]
+    --
+    -- >>> konst i (3::Int,4::Int)
+    -- (3><4)
+    --  [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
+    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
+    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
+    --
+    konst :: e -> d -> c e
+
+instance Container Vector e => Konst e Int Vector
+  where
+    konst = konst'
+
+instance Container Vector e => Konst e (Int,Int) Matrix
+  where
+    konst = konst'
+
+--------------------------------------------------------------------------------
+
+class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
+  where
+    -- |
+    -- >>> build 5 (**2) :: Vector Double
+    -- fromList [0.0,1.0,4.0,9.0,16.0]
+    --
+    -- Hilbert matrix of order N:
+    --
+    -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
+    -- >>> putStr . dispf 2 $ hilb 3
+    -- 3x3
+    -- 1.00  0.50  0.33
+    -- 0.50  0.33  0.25
+    -- 0.33  0.25  0.20
+    --
+    build :: d -> f -> c e
+
+instance Container Vector e => Build Int (e -> e) Vector e
+  where
+    build = build'
+
+instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
+  where
+    build = build'
+
+--------------------------------------------------------------------------------
+
+-- | alternative unicode symbol (25c7) for 'contraction'
+(◇) :: Contraction a b c => a -> b -> c
+infixl 7 ◇
+(◇) = contraction
+
+-- | dot product: @cdot u v = 'udot' ('conj' u) v@
+dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
+dot u v = udot (conj u) v
+
+--------------------------------------------------------------------------------
+
+optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
+optimiseMult = mconcat
+
+--------------------------------------------------------------------------------
+
+
+{- | Compute mean vector and covariance matrix of the rows of a matrix.
+
+>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
+(fromList [4.010341078059521,5.0197204699640405],
+(2><2)
+ [     1.9862461923890056, -1.0127225830525157e-2
+ , -1.0127225830525157e-2,     3.0373954915729318 ])
+
+-}
+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)
+
+--------------------------------------------------------------------------------
+
diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs
index 53fe563..31dcf47 100644
--- a/packages/base/src/Data/Packed/Vector.hs
+++ b/packages/base/src/Data/Packed/Vector.hs
@@ -17,17 +17,15 @@
 
 module Data.Packed.Vector (
     Vector,
-    fromList, (|>), toList, buildVector, constant,
+    fromList, (|>), toList, buildVector,
     dim, (@>),
     subVector, takesV, vjoin, join,
     mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith,
     mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
-    foldLoop, foldVector, foldVectorG, foldVectorWithIndex,
-    stepD, stepF, condD, condF, conjugateC, conjugateQ
+    foldLoop, foldVector, foldVectorG, foldVectorWithIndex
 ) where
 
 import Data.Packed.Internal.Vector
-import Data.Packed.Internal.Matrix
 import Foreign.Storable
 
 -------------------------------------------------------------------
@@ -95,13 +93,3 @@ unzipVector = unzipVectorWith id
 join ::  Storable t => [Vector t] -> Vector t
 join = vjoin
 
-{- | creates a vector with a given number of equal components:
-
-@> constant 2 7
-7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
--}
-constant :: Element a => a -> Int -> Vector a
--- constant x n = runSTVector (newVector x n)
-constant = constantD-- about 2x faster
-
-
diff --git a/packages/base/src/Numeric/Container.hs b/packages/base/src/Numeric/Container.hs
index 067c5fa..6a841aa 100644
--- a/packages/base/src/Numeric/Container.hs
+++ b/packages/base/src/Numeric/Container.hs
@@ -1,258 +1,49 @@
-{-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FunctionalDependencies #-}
-{-# LANGUAGE UndecidableInstances #-}
-
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.Container
--- Copyright   :  (c) Alberto Ruiz 2010-14
--- License     :  BSD3
--- Maintainer  :  Alberto Ruiz
--- Stability   :  provisional
---
--- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
---
--- The 'Container' class is used to define optimized generic functions which work
--- on 'Vector' and 'Matrix' with real or complex elements.
---
--- Some of these functions are also available in the instances of the standard
--- numeric Haskell classes provided by "Numeric.LinearAlgebra".
---
------------------------------------------------------------------------------
 {-# OPTIONS_HADDOCK hide #-}
 
-module Numeric.Container (
-    -- * Basic functions
+module Numeric.Container(
     module Data.Packed,
-    konst, build,
+    constant,
     linspace,
-    diag, ident,
+    diag,
+    ident,
     ctrans,
-    -- * Generic operations
-    Container,
-    add, mul, sub, divide, equal, scaleRecip, addConstant,
-    scalar, conj, scale, arctan2, cmap,
-    atIndex, minIndex, maxIndex, minElement, maxElement,
+    Container(scaleRecip, addConstant,add, sub, mul, divide, equal),
+    scalar,
+    conj,
+    scale,
+    arctan2,
+    cmap,
+    Konst(..),
+    Build(..),
+    atIndex,
+    minIndex, maxIndex, minElement, maxElement,
     sumElements, prodElements,
     step, cond, find, assoc, accum,
-    Transposable(..), Linear(..),
-    -- * Matrix product
-    Product(..), udot, dot, (◇),
-    Mul(..),
-    Contraction(..),(<.>),
+    Element(..),
+    Product(..),
     optimiseMult,
-    mXm,mXv,vXm,LSDiv(..),
+    mXm, mXv, vXm, (<.>),
+    Mul(..),
+    LSDiv, (<\>),
     outer, kronecker,
-    -- * Random numbers
     RandDist(..),
-    randomVector,
-    gaussianSample,
-    uniformSample,
-    -- * Element conversion
+    randomVector, gaussianSample, uniformSample,
+    meanCov,
     Convert(..),
-    Complexable(),
-    RealElement(),
-
-    RealOf, ComplexOf, SingleOf, DoubleOf,
-
-    IndexOf,
+    Complexable,
+    RealElement,
+    RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf,
     module Data.Complex,
-    -- * IO
-    module Data.Packed.IO,
-    -- * Misc
-    Testable(..)
+    dispf, disps, dispcf, vecdisp, latexFormat, format,
+    loadMatrix, saveMatrix, readMatrix
 ) where
 
-import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
-import Data.Packed.Internal.Numeric
-import Data.Complex
-import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
-import Data.Monoid(Monoid(mconcat))
-import Data.Packed.IO
-import Numeric.LinearAlgebra.Random
-
-------------------------------------------------------------------
-
-{- | Creates a real vector containing a range of values:
-
->>> linspace 5 (-3,7::Double)
-fromList [-3.0,-0.5,2.0,4.5,7.0]@
-
->>> linspace 5 (8,2+i) :: Vector (Complex Double)
-fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
-
-Logarithmic spacing can be defined as follows:
-
-@logspace n (a,b) = 10 ** linspace n (a,b)@
--}
-linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
-linspace 0 (a,b) = fromList[(a+b)/2]
-linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
-    where s = (b-a)/fromIntegral (n-1)
-
---------------------------------------------------------
-
-{- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction')
-
-(This operator can also be written using the unicode symbol ◇ (25c7).)
-
-Examples:
-
->>> let a = (3><4)   [1..]      :: Matrix Double
->>> let v = fromList [1,0,2,-1] :: Vector Double
->>> let u = fromList [1,2,3]    :: Vector Double
-
->>> a
-(3><4)
- [ 1.0,  2.0,  3.0,  4.0
- , 5.0,  6.0,  7.0,  8.0
- , 9.0, 10.0, 11.0, 12.0 ]
-
-matrix × matrix:
-
->>> disp 2 (a <.> trans a)
-3x3
- 30   70  110
- 70  174  278
-110  278  446
-
-matrix × vector:
-
->>> a <.> v
-fromList [3.0,11.0,19.0]
-
-dot product:
-
->>> u <.> fromList[3,2,1::Double]
-10
-
-For complex vectors the first argument is conjugated:
-
->>> fromList [1,i] <.> fromList[2*i+1,3]
-1.0 :+ (-1.0)
-
->>> fromList [1,i,1-i] <.> complex a
-fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
--}
-infixl 7 <.>
-(<.>) :: Contraction a b c => a -> b -> c
-(<.>) = contraction
-
 
-class Contraction a b c | a b -> c
-  where
-    -- | Matrix product, matrix - vector product, and dot product
-    contraction :: a -> b -> c
-
-instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
-    u `contraction` v = conj u `udot` v
-
-instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
-    contraction = mXv
-
-instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
-    contraction v m = (conj v) `vXm` m
-
-instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
-    contraction = mXm
-
-
---------------------------------------------------------------------------------
-
-class Mul a b c | a b -> c where
- infixl 7 <>
- -- | Matrix-matrix, matrix-vector, and vector-matrix products.
- (<>)  :: Product t => a t -> b t -> c t
-
-instance Mul Matrix Matrix Matrix where
-    (<>) = mXm
-
-instance Mul Matrix Vector Vector where
-    (<>) m v = flatten $ m <> asColumn v
-
-instance Mul Vector Matrix Vector where
-    (<>) v m = flatten $ asRow v <> m
-
---------------------------------------------------------------------------------
-
-class LSDiv c where
- infixl 7 <\>
- -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
- (<\>)  :: Field t => Matrix t -> c t -> c t
-
-instance LSDiv Vector where
-    m <\> v = flatten (linearSolveSVD m (reshape 1 v))
-
-instance LSDiv Matrix where
-    (<\>) = linearSolveSVD
-
---------------------------------------------------------------------------------
-
-class Konst e d c | d -> c, c -> d
-  where
-    -- |
-    -- >>> konst 7 3 :: Vector Float
-    -- fromList [7.0,7.0,7.0]
-    --
-    -- >>> konst i (3::Int,4::Int)
-    -- (3><4)
-    --  [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
-    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
-    --  , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
-    --
-    konst :: e -> d -> c e
-
-instance Container Vector e => Konst e Int Vector
-  where
-    konst = konst'
-
-instance Container Vector e => Konst e (Int,Int) Matrix
-  where
-    konst = konst'
-
---------------------------------------------------------------------------------
-
-class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
-  where
-    -- |
-    -- >>> build 5 (**2) :: Vector Double
-    -- fromList [0.0,1.0,4.0,9.0,16.0]
-    --
-    -- Hilbert matrix of order N:
-    --
-    -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
-    -- >>> putStr . dispf 2 $ hilb 3
-    -- 3x3
-    -- 1.00  0.50  0.33
-    -- 0.50  0.33  0.25
-    -- 0.33  0.25  0.20
-    --
-    build :: d -> f -> c e
-
-instance Container Vector e => Build Int (e -> e) Vector e
-  where
-    build = build'
-
-instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
-  where
-    build = build'
-
---------------------------------------------------------------------------------
-
--- | alternative unicode symbol (25c7) for 'contraction'
-(◇) :: Contraction a b c => a -> b -> c
-infixl 7 ◇
-(◇) = contraction
-
--- | dot product: @cdot u v = 'udot' ('conj' u) v@
-dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
-dot u v = udot (conj u) v
-
---------------------------------------------------------------------------------
+import Data.Packed.Numeric
+import Data.Packed
+import Data.Packed.Internal(constantD)
+import Data.Complex
 
-optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
-optimiseMult = mconcat
+constant :: Element a => a -> Int -> Vector a
+constant = constantD
 
diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
new file mode 100644
index 0000000..a56c3d2
--- /dev/null
+++ b/packages/base/src/Numeric/HMatrix.hs
@@ -0,0 +1,158 @@
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.HMatrix
+Copyright   :  (c) Alberto Ruiz 2006-14
+License     :  BSD3
+Maintainer  :  Alberto Ruiz
+Stability   :  provisional
+
+-}
+-----------------------------------------------------------------------------
+module Numeric.HMatrix (
+
+    -- * Basic types and data processing    
+    module Numeric.LinearAlgebra.Data,
+    
+    -- * Arithmetic and numeric classes
+    -- |
+    -- The standard numeric classes are defined elementwise:
+    --
+    -- >>> fromList [1,2,3] * fromList [3,0,-2 :: Double]
+    -- fromList [3.0,0.0,-6.0]
+    -- 
+    -- >>> (3><3) [1..9] * ident 3 :: Matrix Double
+    -- (3><3)
+    --  [ 1.0, 0.0, 0.0
+    --  , 0.0, 5.0, 0.0
+    --  , 0.0, 0.0, 9.0 ]
+    --
+    -- In arithmetic operations single-element vectors and matrices
+    -- (created from numeric literals or using 'scalar') automatically
+    -- expand to match the dimensions of the other operand:
+    -- 
+    -- >>> 5 + 2*ident 3 :: Matrix Double
+    -- (3><3)
+    --  [ 7.0, 5.0, 5.0
+    --  , 5.0, 7.0, 5.0
+    --  , 5.0, 5.0, 7.0 ]
+    --
+
+    -- * Matrix product
+    (<.>),
+
+    -- | The overloaded multiplication operators may need type annotations to remove
+    -- ambiguity. In those cases we can also use the specific functions 'mXm', 'mXv', and 'dot'.
+    --
+    -- The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
+    -- single-element matrices (created from numeric literals or using 'scalar')
+    -- are used for scaling.
+    --
+    -- >>> let m = (2><3)[1..] :: Matrix Double
+    -- >>> m <> 2 <> diagl[0.5,1,0]
+    -- (2><3)
+    -- [ 1.0,  4.0, 0.0
+    -- , 4.0, 10.0, 0.0 ]
+    --
+    -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
+
+     
+    -- * Other products
+    outer, kronecker, cross,
+    scale,
+    sumElements, prodElements,
+    
+    -- * Linear Systems
+    (<\>),
+    linearSolve,
+    linearSolveLS,
+    linearSolveSVD,
+    luSolve,
+    cholSolve,
+    cgSolve,
+    cgSolve',
+    
+    -- * Inverse and pseudoinverse
+    inv, pinv, pinvTol,
+
+    -- * Determinant and rank
+    rcond, rank, ranksv, 
+    det, invlndet,
+    
+    -- * Singular value decomposition
+    svd,
+    fullSVD,
+    thinSVD,
+    compactSVD,
+    singularValues,
+    leftSV, rightSV,
+    
+    -- * Eigensystems
+    eig, eigSH, eigSH',
+    eigenvalues, eigenvaluesSH, eigenvaluesSH',
+    geigSH',
+
+    -- * QR
+    qr, rq, qrRaw, qrgr,
+
+    -- * Cholesky
+    chol, cholSH, mbCholSH,
+
+    -- * Hessenberg
+    hess,
+
+    -- * Schur
+    schur,
+
+    -- * LU
+    lu, luPacked,
+    
+    -- * Matrix functions
+    expm,
+    sqrtm,
+    matFunc,
+
+    -- * Nullspace
+    nullspacePrec,
+    nullVector,
+    nullspaceSVD,
+    null1, null1sym,
+    
+    orth,
+
+    -- * Norms
+    norm_0, norm_1, norm_2, norm_Inf,
+    mnorm_0, mnorm_1, mnorm_2, mnorm_Inf,
+    norm_Frob, norm_nuclear,
+
+    -- * Correlation and convolution
+    corr, conv, corrMin, corr2, conv2,
+
+    -- * Random arrays
+
+    RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
+    
+    -- * Misc
+    meanCov, peps, relativeError, haussholder, optimiseMult, dot, udot, mXm, mXv, smXv, (<>), (◇), Seed, checkT,
+    -- * Auxiliary classes
+    Element, Container, Product, Contraction, LSDiv,
+    Complexable, RealElement,
+    RealOf, ComplexOf, SingleOf, DoubleOf,
+    IndexOf,
+    Field,
+    Normed,
+    CGMat, Transposable,
+    ℕ,ℤ,ℝ,ℂ,ℝn,ℂn, 𝑖, i_C --ℍ
+) where
+
+import Numeric.LinearAlgebra.Data
+
+import Numeric.Matrix()
+import Numeric.Vector()
+import Data.Packed.Numeric
+import Numeric.LinearAlgebra.Algorithms
+import Numeric.LinearAlgebra.Util
+import Numeric.LinearAlgebra.Random
+import Numeric.Sparse(smXv)
+import Numeric.LinearAlgebra.Util.CG
+
+
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 242122f..ad315e4 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -1,4 +1,4 @@
------------------------------------------------------------------------------
+--------------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra
 Copyright   :  (c) Alberto Ruiz 2006-14
@@ -7,149 +7,16 @@ Maintainer  :  Alberto Ruiz
 Stability   :  provisional
 
 -}
------------------------------------------------------------------------------
-module Numeric.LinearAlgebra (
-
-    -- * Basic types and data processing    
-    module Numeric.LinearAlgebra.Data,
-    
-    -- * Arithmetic and numeric classes
-    -- |
-    -- The standard numeric classes are defined elementwise:
-    --
-    -- >>> fromList [1,2,3] * fromList [3,0,-2 :: Double]
-    -- fromList [3.0,0.0,-6.0]
-    -- 
-    -- >>> (3><3) [1..9] * ident 3 :: Matrix Double
-    -- (3><3)
-    --  [ 1.0, 0.0, 0.0
-    --  , 0.0, 5.0, 0.0
-    --  , 0.0, 0.0, 9.0 ]
-    --
-    -- In arithmetic operations single-element vectors and matrices
-    -- (created from numeric literals or using 'scalar') automatically
-    -- expand to match the dimensions of the other operand:
-    -- 
-    -- >>> 5 + 2*ident 3 :: Matrix Double
-    -- (3><3)
-    --  [ 7.0, 5.0, 5.0
-    --  , 5.0, 7.0, 5.0
-    --  , 5.0, 5.0, 7.0 ]
-    --
-
-    -- * Matrix product
-    (<.>),
-
-    -- | The overloaded multiplication operators may need type annotations to remove
-    -- ambiguity. In those cases we can also use the specific functions 'mXm', 'mXv', and 'dot'.
-    --
-    -- The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
-    -- single-element matrices (created from numeric literals or using 'scalar')
-    -- are used for scaling.
-    --
-    -- >>> let m = (2><3)[1..] :: Matrix Double
-    -- >>> m <> 2 <> diagl[0.5,1,0]
-    -- (2><3)
-    -- [ 1.0,  4.0, 0.0
-    -- , 4.0, 10.0, 0.0 ]
-    --
-    -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
-
-     
-    -- * Other products
-    outer, kronecker, cross,
-    scale,
-    sumElements, prodElements, absSum,
-    
-    -- * Linear Systems
-    (<\>),
-    linearSolve,
-    linearSolveLS,
-    linearSolveSVD,
-    luSolve,
-    cholSolve,
-    cgSolve,
-    cgSolve',
-    
-    -- * Inverse and pseudoinverse
-    inv, pinv, pinvTol,
-
-    -- * Determinant and rank
-    rcond, rank, ranksv, 
-    det, invlndet,
-    
-    -- * Singular value decomposition
-    svd,
-    fullSVD,
-    thinSVD,
-    compactSVD,
-    singularValues,
-    leftSV, rightSV,
-    
-    -- * Eigensystems
-    eig, eigSH, eigSH',
-    eigenvalues, eigenvaluesSH, eigenvaluesSH',
-    geigSH',
-
-    -- * QR
-    qr, rq, qrRaw, qrgr,
-
-    -- * Cholesky
-    chol, cholSH, mbCholSH,
-
-    -- * Hessenberg
-    hess,
-
-    -- * Schur
-    schur,
-
-    -- * LU
-    lu, luPacked,
-    
-    -- * Matrix functions
-    expm,
-    sqrtm,
-    matFunc,
+--------------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
 
-    -- * Nullspace
-    nullspacePrec,
-    nullVector,
-    nullspaceSVD,
-    null1, null1sym,
-    
-    orth,
-
-    -- * Norms
-    norm1, norm2, normInf, pnorm, NormType(..),
-
-    -- * Correlation and convolution
-    corr, conv, corrMin, corr2, conv2,
-
-    -- * Random arrays
-
-    RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
-    
-    -- * Misc
-    meanCov, peps, relativeError, haussholder, optimiseMult, dot, udot, mXm, mXv, smXv, (<>), (◇), Seed, checkT,
-    -- * Auxiliary classes
-    Element, Container, Product, Contraction, LSDiv,
-    Complexable(), RealElement(),
-    RealOf, ComplexOf, SingleOf, DoubleOf,
-    IndexOf,
-    Field, Normed,
-    CGMat, Transposable,
-    R,V
+module Numeric.LinearAlgebra (
+    module Numeric.Container,
+    module Numeric.LinearAlgebra.Algorithms
 ) where
 
-import Numeric.LinearAlgebra.Data
-
-import Numeric.Matrix()
-import Numeric.Vector()
 import Numeric.Container
 import Numeric.LinearAlgebra.Algorithms
-import Numeric.LinearAlgebra.Util
-import Numeric.LinearAlgebra.Random
-import Numeric.Sparse(smXv)
-import Numeric.LinearAlgebra.Util.CG
-
+import Numeric.Matrix()
+import Numeric.Vector()
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
index c7e7043..bbcc513 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -76,12 +76,12 @@ module Numeric.LinearAlgebra.Algorithms (
 ) where
 
 
-import Data.Packed.Development hiding ((//))
 import Data.Packed
 import Numeric.LinearAlgebra.LAPACK as LAPACK
 import Data.List(foldl1')
 import Data.Array
 import Data.Packed.Internal.Numeric
+import Data.Packed.Internal(shSize)
 
 
 {- | Generic linear algebra functions for double precision real and complex matrices.
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index e3cbe31..89bebbe 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -16,13 +16,19 @@ module Numeric.LinearAlgebra.Data(
     -- * Vector
     -- | 1D arrays are storable vectors from the vector package.
     
-    Vector, (|>), dim, (@>),
-    
+    vect, (|>),
+
     -- * Matrix
-    Matrix, (><), size, (@@>), trans, ctrans,
+    
+    mat, (><), tr,
+    
+    -- * Indexing
+    
+    size,
+    Indexable(..),
     
     -- * Construction
-    scalar, konst, build, assoc, accum, linspace, -- ones, zeros,
+    scalar, Konst(..), Build(..), assoc, accum, linspace, -- ones, zeros,
 
     -- * Diagonal
     ident, diag, diagl, diagRect, takeDiag,
@@ -62,11 +68,13 @@ module Numeric.LinearAlgebra.Data(
 
     module Data.Complex,
 
+    Vector, Matrix
+
 ) where
 
 import Data.Packed.Vector
 import Data.Packed.Matrix
-import Numeric.Container
+import Data.Packed.Numeric
 import Numeric.LinearAlgebra.Util
 import Data.Complex
 import Numeric.Sparse
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
index 2f91e18..a7d6946 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -1,4 +1,9 @@
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Util
@@ -14,19 +19,24 @@ Stability   :  provisional
 module Numeric.LinearAlgebra.Util(
 
     -- * Convenience functions
-    size, disp,
+    vect, mat,
+    disp,
     zeros, ones,
     diagl,
     row,
     col,
     (&), (¦), (——), (#),
     (?), (¿),
+    Indexable(..), size,
+    rand, randn,
     cross,
     norm,
+    ℕ,ℤ,ℝ,ℂ,ℝn,ℂn,𝑖,i_C, --ℍ
+    norm_1, norm_2, norm_0, norm_Inf, norm_Frob, norm_nuclear,
+    mnorm_1, mnorm_2, mnorm_0, mnorm_Inf,
     unitary,
     mt,
     pairwiseD2,
-    meanCov,
     rowOuters,
     null1,
     null1sym,
@@ -48,13 +58,49 @@ module Numeric.LinearAlgebra.Util(
     vtrans
 ) where
 
-import Numeric.Container
+import Data.Packed.Numeric
 import Numeric.LinearAlgebra.Algorithms hiding (i)
 import Numeric.Matrix()
 import Numeric.Vector()
-
+import Numeric.LinearAlgebra.Random
 import Numeric.LinearAlgebra.Util.Convolution
 
+type ℝ = Double
+type ℕ = Int
+type ℤ = Int
+type ℂ = Complex Double
+type ℝn = Vector ℝ
+type ℂn = Vector ℂ
+--newtype ℍ m = H m
+
+i_C, 𝑖 :: ℂ
+𝑖 = 0:+1
+i_C = 𝑖
+
+{- | create a real vector
+
+>>> vect [1..5]
+fromList [1.0,2.0,3.0,4.0,5.0]
+
+-}
+vect :: [ℝ] -> ℝn
+vect = fromList
+
+{- | create a real matrix
+
+>>> mat 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 ]
+
+-}
+mat
+  :: Int -- ^ columns
+  -> [ℝ] -- ^ elements
+  -> Matrix ℝ
+mat c = reshape c . fromList
+
 {- | print a real matrix with given number of digits after the decimal point
 
 >>> disp 5 $ ident 2 / 3
@@ -175,38 +221,97 @@ norm :: Vector Double -> Double
 -- ^ 2-norm of real vector
 norm = pnorm PNorm2
 
+norm_2 :: Normed Vector t => Vector t -> RealOf t
+norm_2 = pnorm PNorm2
+
+norm_1 :: Normed Vector t => Vector t -> RealOf t
+norm_1 = pnorm PNorm1
+
+norm_Inf :: Normed Vector t => Vector t -> RealOf t
+norm_Inf = pnorm Infinity
+
+norm_0 :: Vector ℝ -> ℝ
+norm_0 v = sumElements (step (abs v - scalar (eps*mx)))
+  where
+    mx = norm_Inf v
+
+norm_Frob :: Normed Matrix t => Matrix t -> RealOf t
+norm_Frob = pnorm Frobenius
+
+norm_nuclear :: Field t => Matrix t -> ℝ
+norm_nuclear = sumElements . singularValues
+
+mnorm_2 :: Normed Matrix t => Matrix t -> RealOf t
+mnorm_2 = pnorm PNorm2
+
+mnorm_1 :: Normed Matrix t => Matrix t -> RealOf t
+mnorm_1 = pnorm PNorm1
+
+mnorm_Inf :: Normed Matrix t => Matrix t -> RealOf t
+mnorm_Inf = pnorm Infinity
+
+mnorm_0 :: Matrix ℝ -> ℝ
+mnorm_0 = norm_0 . flatten
 
 -- | Obtains a vector in the same direction with 2-norm=1
 unitary :: Vector Double -> Vector Double
 unitary v = v / scalar (norm v)
 
--- | ('rows' &&& 'cols')
-size :: Matrix t -> (Int, Int)
-size m = (rows m, cols m)
 
 -- | trans . inv
 mt :: Matrix Double -> Matrix Double
 mt = trans . inv
 
 --------------------------------------------------------------------------------
+{- |
+
+>>> size $ fromList[1..10::Double]
+10
+>>> size $ (2><5)[1..10::Double]
+(2,5)
+
+-}
+size :: Container c t => c t -> IndexOf c
+size = size'
 
-{- | Compute mean vector and covariance matrix of the rows of a matrix.
+{- |
+
+>>> vect [1..10] ! 3
+4.0
+
+>>> mat 5 [1..15] ! 1
+fromList [6.0,7.0,8.0,9.0,10.0]
 
->>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
-(fromList [4.010341078059521,5.0197204699640405],
-(2><2)
- [     1.9862461923890056, -1.0127225830525157e-2
- , -1.0127225830525157e-2,     3.0373954915729318 ])
+>>> mat 5 [1..15] ! 1 ! 3
+9.0
 
 -}
-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)
+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 (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
 
 --------------------------------------------------------------------------------
 
@@ -220,7 +325,7 @@ pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
     ox = one (rows x)
     oy = one (rows y)
     oc = one (cols x)
-    one k = constant 1 k
+    one k = konst 1 k
     x2 = x * x <> oc
     y2 = y * y <> oc
     ok = cols x == cols y
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
index d21602d..5e2ea84 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
@@ -6,7 +6,7 @@ module Numeric.LinearAlgebra.Util.CG(
     CGMat, CGState(..), R, V
 ) where
 
-import Numeric.Container
+import Data.Packed.Numeric
 import Numeric.Vector()
 
 {-
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
index e4cba8f..c8c7536 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -16,16 +16,18 @@ module Numeric.LinearAlgebra.Util.Convolution(
    corr2, conv2, separable
 ) where
 
-import Numeric.Container
+import Data.Packed.Numeric
 
 
 vectSS :: Element t => Int -> Vector t -> Matrix t
 vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ]
 
 
-corr :: Product t => Vector t -- ^ kernel
-                  -> Vector t -- ^ source
-                  -> Vector t
+corr
+  :: (Container Vector t, Product t)
+    => Vector t -- ^ kernel
+    -> Vector t -- ^ source
+    -> Vector t
 {- ^ correlation
 
 >>> corr (fromList[1,2,3]) (fromList [1..10])
@@ -33,12 +35,12 @@ fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]
 
 -}
 corr ker v
-    | dim ker == 0 = constant 0 (dim v)
+    | dim ker == 0 = konst 0 (dim v)
     | dim ker <= dim v = vectSS (dim ker) v <> ker
     | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")"
 
 
-conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t
+conv :: (Container Vector t, Product t, Num t) => Vector t -> Vector t -> Vector t
 {- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product)
 
 >>> conv (fromList[1,1]) (fromList [-1,1])
@@ -46,12 +48,12 @@ fromList [-1.0,0.0,1.0]
 
 -}
 conv ker v
-    | dim ker == 0 = constant 0 (dim v)
+    | dim ker == 0 = konst 0 (dim v)
     | otherwise = corr ker' v'
   where
     ker' = (flatten.fliprl.asRow) ker
     v' = vjoin [z,v,z]
-    z = constant 0 (dim ker -1)
+    z = konst 0 (dim ker -1)
 
 corrMin :: (Container Vector t, RealElement t, Product t)
         => Vector t
diff --git a/packages/base/src/Numeric/Sparse.hs b/packages/base/src/Numeric/Sparse.hs
index 1957d3a..2df4578 100644
--- a/packages/base/src/Numeric/Sparse.hs
+++ b/packages/base/src/Numeric/Sparse.hs
@@ -10,7 +10,7 @@ module Numeric.Sparse(
     smXv
 )where
 
-import Numeric.Container
+import Data.Packed.Numeric
 import qualified Data.Vector.Storable as V
 import Data.Function(on)
 import Control.Arrow((***))
@@ -133,10 +133,6 @@ instance Transposable SMatrix
     tr (CSC vs rs cs n m) = CSR vs rs cs m n
     tr (Diag v n m) = Diag v m n
 
-instance Transposable (Matrix Double)
-  where
-    tr = trans
-
 
 instance CGMat SMatrix
 instance CGMat (Matrix Double)