extend conformability to empty arrays

3 files changed, 69 insertions, 19 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 []