10 files changed, 133 insertions, 59 deletions

diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index 49f17b0..edef3c2 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -49,7 +49,11 @@ common f = commonval . map f where
 compatdim :: [Int] -> Maybe Int
 compatdim [] = Nothing
 compatdim [a] = Just a
-compatdim (a:b:xs) = if a==b || a==1 || b==1 then compatdim (max a b:xs) else Nothing
+compatdim (a:b:xs)
+    | a==b = compatdim (b:xs)
+    | a==1 = compatdim (b:xs)
+    | b==1 = compatdim (a:xs)
+    | otherwise = Nothing
 
 -- | Formatting tool
 table :: String -> [[String]] -> String
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 2004e85..9719fc0 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -32,6 +32,7 @@ module Data.Packed.Internal.Matrix(
     (@@>), atM',
     saveMatrix,
     singleton,
+    emptyM,
     size, shSize, conformVs, conformMs, conformVTo, conformMTo
 ) where
 
@@ -157,16 +158,24 @@ toLists m = splitEvery (cols m) . toList . flatten $ m
 -- All vectors must have the same dimension,
 -- or dimension 1, which is are automatically expanded.
 fromRows :: Element t => [Vector t] -> Matrix t
+fromRows [] = emptyM 0 0
 fromRows vs = case compatdim (map dim vs) of
-    Nothing -> error "fromRows applied to [] or to vectors with different sizes"
-    Just c  -> reshape c . vjoin . map (adapt c) $ vs
+    Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs)
+    Just 0  -> emptyM r 0
+    Just c  -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs
   where
-    adapt c v | dim v == c = v
-              | otherwise = constantD (v@>0) c
+    r = length vs
+    adapt c v
+        | c == 0 = fromList[]
+        | dim v == c = v
+        | otherwise = constantD (v@>0) c
 
 -- | extracts the rows of a matrix as a list of vectors
 toRows :: Element t => Matrix t -> [Vector t]
-toRows m = toRows' 0 where
+toRows m
+    | c == 0    = replicate r (fromList[])
+    | otherwise = toRows' 0
+  where
     v = flatten m
     r = rows m
     c = cols m
@@ -200,7 +209,7 @@ atM' Matrix {irows = r, xdat = v, order = ColumnMajor} i j = v `at'` (j*r+i)
 
 matrixFromVector o r c v
     | r * c == dim v = m
-    | otherwise = error $ "matrixFromVector " ++ shSize m ++ " <- " ++ show (dim v)
+    | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
   where
     m = Matrix { irows = r, icols = c, xdat = v, order = o }
 
@@ -398,8 +407,8 @@ subMatrix :: Element a
           -> Matrix a -- ^ input matrix
           -> Matrix a -- ^ result
 subMatrix (r0,c0) (rt,ct) m
-    | 0 <= r0 && 0 < rt && r0+rt <= (rows m) &&
-      0 <= c0 && 0 < ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
+    | 0 <= r0 && 0 <= rt && r0+rt <= (rows m) &&
+      0 <= c0 && 0 <= ct && c0+ct <= (cols m) = subMatrixD (r0,c0) (rt,ct) m
     | otherwise = error $ "wrong subMatrix "++
                           show ((r0,c0),(rt,ct))++" of "++show(rows m)++"x"++ show (cols m)
 
@@ -437,18 +446,21 @@ foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr
 
 ----------------------------------------------------------------------
 
+maxZ xs = if minimum xs == 0 then 0 else maximum xs
+
 conformMs ms = map (conformMTo (r,c)) ms
   where
-    r = maximum (map rows ms)
-    c = maximum (map cols ms)
+    r = maxZ (map rows ms)
+    c = maxZ (map cols ms)
+    
 
 conformVs vs = map (conformVTo n) vs
   where
-    n = maximum (map dim vs)
+    n = maxZ (map dim vs)
 
 conformMTo (r,c) m
     | size m == (r,c) = m
-    | size m == (1,1) = reshape c (constantD (m@@>(0,0)) (r*c))
+    | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
     | size m == (r,1) = repCols c m
     | size m == (1,c) = repRows r m
     | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to (" ++ show r ++ "><"++ show c ++")"
@@ -465,6 +477,8 @@ size m = (rows m, cols m)
 
 shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
 
+emptyM r c = matrixFromVector RowMajor r c (fromList[])
+
 ----------------------------------------------------------------------
 
 instance (Storable t, NFData t) => NFData (Matrix t)
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index b92d60f..d94d167 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -35,7 +35,7 @@ module Data.Packed.Matrix (
     repmat,
     flipud, fliprl,
     subMatrix, takeRows, dropRows, takeColumns, dropColumns,
-    extractRows,
+    extractRows, extractColumns,
     diagRect, takeDiag,
     mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
     liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D
@@ -104,6 +104,7 @@ breakAt c l = (a++[c],tail b) where
 
 -- | creates a matrix from a vertical list of matrices
 joinVert :: Element t => [Matrix t] -> Matrix t
+joinVert [] = emptyM 0 0
 joinVert ms = case common cols ms of
     Nothing -> error "(impossible) joinVert on matrices with different number of columns"
     Just c  -> matrixFromVector RowMajor (sum (map rows ms)) c $ vjoin (map flatten ms)
@@ -173,6 +174,11 @@ adaptBlocks ms = ms' where
 0  0  0  0  0  0  0  5
 0  0  0  0  0  0  0  7
 
+>>> diagBlock [(0><4)[], konst 2 (2,3)]  :: Matrix Double
+(2><7)
+ [ 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0
+ , 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0 ]
+
 -}
 diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
 diagBlock ms = fromBlocks $ zipWith f ms [0..]
@@ -186,11 +192,15 @@ diagBlock ms = fromBlocks $ zipWith f ms [0..]
 
 -- | Reverse rows
 flipud :: Element t => Matrix t -> Matrix t
-flipud m = fromRows . reverse . toRows $ m
+flipud m = extractRows [r-1,r-2 .. 0] $ m
+  where
+    r = rows m
 
 -- | Reverse columns
 fliprl :: Element t => Matrix t -> Matrix t
-fliprl m = fromColumns . reverse . toColumns $ m
+fliprl m = extractColumns [c-1,c-2 .. 0] $ m
+  where
+    c = cols m
 
 ------------------------------------------------------------
 
@@ -327,8 +337,25 @@ fromArray2D m = (r><c) (elems m)
 
 -- | rearranges the rows of a matrix according to the order given in a list of integers.
 extractRows :: Element t => [Int] -> Matrix t -> Matrix t
+extractRows [] m = emptyM 0 (cols m)
 extractRows l m = fromRows $ extract (toRows m) l
-    where extract l' is = [l'!!i |i<-is]
+  where
+    extract l' is = [l'!!i | i<- map verify is]
+    verify k
+        | k >= 0 && k < rows m = k
+        | otherwise = error $ "can't extract row "
+                           ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
+
+-- | rearranges the rows of a matrix according to the order given in a list of integers.
+extractColumns :: Element t => [Int] -> Matrix t -> Matrix t
+extractColumns l m = trans . extractRows (map verify l) . trans $ m
+  where
+    verify k
+        | k >= 0 && k < cols m = k
+        | otherwise = error $ "can't extract column "
+                           ++show k++" in list " ++ show l ++ " from matrix " ++ shSize m
+
+
 
 {- | creates matrix by repetition of a matrix a given number of rows and columns
 
@@ -341,7 +368,9 @@ extractRows l m = fromRows $ extract (toRows m) l
 
 -}
 repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
-repmat m r c = fromBlocks $ splitEvery c $ replicate (r*c) m
+repmat m r c
+    | r == 0 || c == 0 = emptyM (r*rows m) (c*cols m)
+    | otherwise = fromBlocks $ replicate r $ replicate c $ m
 
 -- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.
 liftMatrix2Auto :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
@@ -390,7 +419,10 @@ toBlocks rs cs m = map (toBlockCols cs) . toBlockRows rs $ m
 -- | Fully partition a matrix into blocks of the same size. If the dimensions are not
 -- a multiple of the given size the last blocks will be smaller.
 toBlocksEvery :: (Element t) => Int -> Int -> Matrix t -> [[Matrix t]]
-toBlocksEvery r c m = toBlocks rs cs m where
+toBlocksEvery r c m
+    | r < 1 || c < 1 = error $ "toBlocksEvery expects block sizes > 0, given "++show r++" and "++ show c
+    | otherwise = toBlocks rs cs m
+  where
     (qr,rr) = rows m `divMod` r
     (qc,rc) = cols m `divMod` c
     rs = replicate qr r ++ if rr > 0 then [rr] else []
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index b145a26..dea8a79 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -86,15 +86,19 @@ constant = constantD-- about 2x faster
 
 {- | Creates a real vector containing a range of values:
 
->>> linspace 5 (-3,7)
+>>> 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 :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e
-linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
+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)
 
 -- | dot product: @cdot u v = 'udot' ('conj' u) v@
diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs
index 6445e04..ea4262c 100644
--- a/lib/Numeric/ContainerBoot.hs
+++ b/lib/Numeric/ContainerBoot.hs
@@ -45,7 +45,7 @@ import Numeric.Conversion
 import Data.Packed.Internal
 import Numeric.GSL.Vector
 import Data.Complex
-import Control.Monad(ap)
+import Control.Applicative((<*>))
 
 import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
 
@@ -206,10 +206,10 @@ instance Container Vector Float where
     conj = id
     cmap = mapVector
     atIndex = (@>)
-    minIndex     = round . toScalarF MinIdx
-    maxIndex     = round . toScalarF MaxIdx
-    minElement  = toScalarF Min
-    maxElement  = toScalarF Max
+    minIndex     = emptyErrorV "minIndex"   (round . toScalarF MinIdx)
+    maxIndex     = emptyErrorV "maxIndex"   (round . toScalarF MaxIdx)
+    minElement   = emptyErrorV "minElement" (toScalarF Min)
+    maxElement   = emptyErrorV "maxElement" (toScalarF Max)
     sumElements  = sumF
     prodElements = prodF
     step = stepF
@@ -234,10 +234,10 @@ instance Container Vector Double where
     conj = id
     cmap = mapVector
     atIndex = (@>)
-    minIndex     = round . toScalarR MinIdx
-    maxIndex     = round . toScalarR MaxIdx
-    minElement  = toScalarR Min
-    maxElement  = toScalarR Max
+    minIndex     = emptyErrorV "minIndex"   (round . toScalarR MinIdx)
+    maxIndex     = emptyErrorV "maxIndex"   (round . toScalarR MaxIdx)
+    minElement   = emptyErrorV "minElement" (toScalarR Min)
+    maxElement   = emptyErrorV "maxElement" (toScalarR Max)
     sumElements  = sumR
     prodElements = prodR
     step = stepD
@@ -262,10 +262,10 @@ instance Container Vector (Complex Double) where
     conj = conjugateC
     cmap = mapVector
     atIndex = (@>)
-    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    minElement  = ap (@>) minIndex
-    maxElement  = ap (@>) maxIndex
+    minIndex     = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj))
+    maxIndex     = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj))
+    minElement   = emptyErrorV "minElement" (atIndex <*> minIndex)
+    maxElement   = emptyErrorV "maxElement" (atIndex <*> maxIndex)
     sumElements  = sumC
     prodElements = prodC
     step = undefined -- cannot match
@@ -290,10 +290,10 @@ instance Container Vector (Complex Float) where
     conj = conjugateQ
     cmap = mapVector
     atIndex = (@>)
-    minIndex     = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    maxIndex     = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
-    minElement  = ap (@>) minIndex
-    maxElement  = ap (@>) maxIndex
+    minIndex     = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj))
+    maxIndex     = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj))
+    minElement   = emptyErrorV "minElement" (atIndex <*> minIndex)
+    maxElement   = emptyErrorV "maxElement" (atIndex <*> maxIndex)
     sumElements  = sumQ
     prodElements = prodQ
     step = undefined -- cannot match
@@ -320,14 +320,12 @@ instance (Container Vector a) => Container Matrix a where
     conj = liftMatrix conj
     cmap f = liftMatrix (mapVector f)
     atIndex = (@@>)
-    minIndex m = let (r,c) = (rows m,cols m)
-                     i = (minIndex $ flatten m)
-                 in (i `div` c,i `mod` c)
-    maxIndex m = let (r,c) = (rows m,cols m)
-                     i = (maxIndex $ flatten m)
-                 in (i `div` c,i `mod` c)
-    minElement = ap (@@>) minIndex
-    maxElement = ap (@@>) maxIndex
+    minIndex = emptyErrorM "minIndex of Matrix" $
+                \m -> divMod (minIndex $ flatten m) (cols m)
+    maxIndex = emptyErrorM "maxIndex of Matrix" $
+                \m -> divMod (maxIndex $ flatten m) (cols m)
+    minElement = emptyErrorM "minElement of Matrix" (atIndex <*> minIndex)
+    maxElement = emptyErrorM "maxElement of Matrix" (atIndex <*> maxIndex)
     sumElements = sumElements . flatten
     prodElements = prodElements . flatten
     step = liftMatrix step
@@ -336,6 +334,17 @@ instance (Container Vector a) => Container Matrix a where
     accum = accumM
     cond = condM
 
+
+emptyErrorV msg f v =
+    if dim v > 0
+        then f v
+        else error $ msg ++ " of Vector with dim = 0"
+
+emptyErrorM msg f m =
+    if rows m > 0 && cols m > 0
+        then f m
+        else error $ msg++" "++shSize m
+
 ----------------------------------------------------
 
 -- | Matrix product and related functions
@@ -393,7 +402,6 @@ emptyVal f v =
         then f v
         else 0
 
-
 -- FIXME remove unused C wrappers
 -- | (unconjugated) dot product
 udot :: Product e => Vector e -> Vector e -> e
@@ -592,7 +600,7 @@ accumM m0 f xs = ST.runSTMatrix $ do
 
 ----------------------------------------------------------------------
 
-condM a b l e t = reshape (cols a'') $ cond a' b' l' e' t'
+condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t'
   where
     args@(a'':_) = conformMs [a,b,l,e,t]
     [a', b', l', e', t'] = map flatten args
diff --git a/lib/Numeric/HMatrix.hs b/lib/Numeric/HMatrix.hs
index a2f09df..f49ea53 100644
--- a/lib/Numeric/HMatrix.hs
+++ b/lib/Numeric/HMatrix.hs
@@ -114,7 +114,7 @@ module Numeric.HMatrix (
     orth,
 
     -- * Norms
-    norm1, norm2, normInf,
+    norm1, norm2, normInf, pnorm, NormType(..),
 
     -- * Correlation and Convolution
     corr, conv, corrMin, corr2, conv2,
diff --git a/lib/Numeric/HMatrix/Data.hs b/lib/Numeric/HMatrix/Data.hs
index 288b0af..568dc05 100644
--- a/lib/Numeric/HMatrix/Data.hs
+++ b/lib/Numeric/HMatrix/Data.hs
@@ -44,7 +44,7 @@ module Numeric.HMatrix.Data(
     find, maxIndex, minIndex, maxElement, minElement, atIndex,
 
     -- * IO
-    disp, dispf, disps, dispcf, vecdisp, latexFormat, format,
+    disp, dispf, disps, dispcf, latexFormat, format,
     loadMatrix, saveMatrix, fromFile, fileDimensions,
     readMatrix,
     fscanfVector, fprintfVector, freadVector, fwriteVector,
diff --git a/lib/Numeric/HMatrix/Devel.hs b/lib/Numeric/HMatrix/Devel.hs
index 7363477..b921f44 100644
--- a/lib/Numeric/HMatrix/Devel.hs
+++ b/lib/Numeric/HMatrix/Devel.hs
@@ -55,7 +55,7 @@ module Numeric.HMatrix.Devel(
     Complexable(), RealElement(),
     RealOf, ComplexOf, SingleOf, DoubleOf,
     IndexOf,
-    Field,
+    Field, Normed
 ) where
 
 import Data.Packed.Foreign
@@ -65,5 +65,5 @@ import Numeric.Container(Container,Contraction,LSDiv,Product,
                          Complexable(),RealElement(),
                          RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf)
 import Data.Packed
-import Numeric.LinearAlgebra.Algorithms(Field)
+import Numeric.LinearAlgebra.Algorithms(Field,Normed)
 
diff --git a/lib/Numeric/IO.hs b/lib/Numeric/IO.hs
index 57275ac..836f352 100644
--- a/lib/Numeric/IO.hs
+++ b/lib/Numeric/IO.hs
@@ -60,6 +60,9 @@ disps d x = sdims x ++ "  " ++ formatScaled d x
 3.00  3.50  4.00  4.50
 5.00  5.50  6.00  6.50
 
+>>> putStr . unlines . tail . lines . dispf 2 . asRow $ linspace 10 (0,1)
+0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
+
 -}
 dispf :: Int -> Matrix Double -> String
 dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
@@ -74,7 +77,8 @@ formatScaled dec t = "E"++show o++"\n" ++ ss
     where ss = format " " (printf fmt. g) t
           g x | o >= 0    = x/10^(o::Int)
               | otherwise = x*10^(-o)
-          o = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
+          o | rows t == 0 || cols t == 0 = 0
+            | otherwise = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
           fmt = '%':show (dec+3) ++ '.':show dec ++"f"
 
 {- | Show a vector using a function for showing matrices.
diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs
index 7164827..7d134bf 100644
--- a/lib/Numeric/LinearAlgebra/Util.hs
+++ b/lib/Numeric/LinearAlgebra/Util.hs
@@ -166,7 +166,7 @@ row = asRow . fromList
 col :: [Double] -> Matrix Double
 col = asColumn . fromList
 
-{- | extract selected rows
+{- | extract rows
 
 >>> (20><4) [1..] ? [2,1,1]
 (3><4)
@@ -179,12 +179,20 @@ infixl 9 ?
 (?) :: Element t => Matrix t -> [Int] -> Matrix t
 (?) = flip extractRows
 
--- | extract selected columns
---
--- (unicode 0x00bf, inverted question mark)
+{- | 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
-m ¿ ks = trans . extractRows ks . trans $ m
+(¿)= flip extractColumns
 
 
 cross :: Vector Double -> Vector Double -> Vector Double