7 files changed, 64 insertions, 40 deletions

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