From c9914d694d3b86ece46fa0c76e0466c6cd394d14 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 6 May 2014 08:50:50 +0200 Subject: extend conformability to empty arrays --- lib/Data/Packed/Internal/Common.hs | 6 +++- lib/Data/Packed/Internal/Matrix.hs | 38 +++++++++++++++-------- lib/Data/Packed/Matrix.hs | 44 +++++++++++++++++++++++---- lib/Numeric/Container.hs | 10 ++++-- lib/Numeric/ContainerBoot.hs | 62 +++++++++++++++++++++----------------- lib/Numeric/HMatrix.hs | 2 +- lib/Numeric/HMatrix/Data.hs | 2 +- lib/Numeric/HMatrix/Devel.hs | 4 +-- lib/Numeric/IO.hs | 6 +++- lib/Numeric/LinearAlgebra/Util.hs | 18 ++++++++--- 10 files changed, 133 insertions(+), 59 deletions(-) (limited to 'lib') 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> [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 -- cgit v1.2.3