dsp 0.2.3.1 → 0.2.4

raw patch · 9 files changed

+759/−49 lines, 9 filesdep +QuickCheckdep +containersdep ~arraydep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: QuickCheck, containers

Dependency ranges changed: array, base

API changes (from Hackage documentation)

Files

• Matrix/LU.hs +35/−3
• Matrix/Matrix.hs +86/−28
• Matrix/QR/Givens.hs +160/−0
• Matrix/QR/Householder.hs +152/−0
• Matrix/Sparse.hs +84/−0
• Matrix/Vector.hs +35/−0
• dsp.cabal +28/−18
• test/Main.hs +15/−0
• test/Test/Matrix/QR.hs +164/−0

@@ -14,8 +14,13 @@  module Matrix.LU (lu, lu_solve, improve, inverse, lu_det, solve, det) where +import qualified Matrix.Matrix as Matrix+import qualified Matrix.Vector as Vector+import qualified Data.List as List import Data.Array+import Data.Ord (comparing) + -- | LU decomposition via Crout's Algorithm  -- TODO: modify for partial pivoting / permutation matrix@@ -37,8 +42,8 @@  -- | Solution to Ax=b via LU decomposition --- forward is forumla (2.3.6) in NRIC, but remebering that a11=1--- backward is forumla (2.3.7) in NRIC+-- forward is formula (2.3.6) in NRIC, but remembering that a11=1+-- backward is formula (2.3.7) in NRIC  lu_solve :: Array (Int,Int) Double -- ^ LU(A)          -> Array Int Double -- ^ b@@ -107,10 +112,37 @@  -- | determinant using original matrix +{-+It is based on LU decomposition without singularity check+and thus returns NaN instead of zero if the matrix is singular.+-}+_det :: Array (Int,Int) Double -- ^ A+    -> Double -- ^ det(A)++_det a = (lu_det . lu) a++{- |+Determinant computation by implicit LU decomposition with permutations.+-} det :: Array (Int,Int) Double -- ^ A     -> Double -- ^ det(A) -det a = (lu_det . lu) a+det a =+   if rangeSize (bounds a) == 0+     then 1+     else+         let ((m0,n0), (m1,n1)) = bounds a+             v = Matrix.getColumn n0 a+             (maxi,maxv) = List.maximumBy (comparing (abs . snd)) $ assocs v+             reduced =+                ixmap ((m0,n0), (pred m1, pred n1))+                   (\(i,j) -> (if i<maxi then i else succ i, succ j)) $+                Vector.sub a $ Matrix.outer v $+                Vector.scale (recip maxv) $ Matrix.getRow maxi a+             sign = if even (rangeSize (m0,maxi)-1) then 1 else -1+             pivot = a!(maxi,n0)+         in  if pivot == 0 then 0 else sign * pivot * det reduced+  ------------------------------------------------------------------------------- -- tests

@@ -14,51 +14,109 @@  module Matrix.Matrix where +import Matrix.Vector (generate) import Data.Array import Data.Complex  -- | Matrix-matrix multiplication: A x B = C -mm_mult :: (Ix a, Integral a, Num b) => Array (a,a) b -- ^ A-	-> Array (a,a) b -- ^ B-	-> Array (a,a) b -- ^ C+mm_mult :: (Ix i, Ix j, Ix k, Num a) => Array (i,j) a -- ^ A+	-> Array (j,k) a -- ^ B+	-> Array (i,k) a -- ^ C -mm_mult a b = if ac /= br+mm_mult a b = if (ac0,ac1) /= (br0,br1) 	      then error "mm_mult: inside dimensions inconsistent"-	      else array bnds [ ((i,j), mult i j) | (i,j) <- range bnds ]-    where mult i j = sum [ a!(i,k) * b!(k,j) | k <- [1..ac] ]-	  ((_,_),(ar,ac)) = bounds a-	  ((_,_),(br,bc)) = bounds b-	  bnds = ((1,1),(ar,bc))+	      else generate ((ar0,bc0),(ar1,bc1)) $ \(i,j) ->+			sum [ a!(i,k) * b!(k,j) | k <- range (ac0,ac1) ]+    where ((ar0,ac0),(ar1,ac1)) = bounds a+	  ((br0,bc0),(br1,bc1)) = bounds b  -- | Matrix-vector multiplication: A x b = c -mv_mult :: (Ix a, Integral a, Num b) => Array (a,a) b -- ^ A-	-> Array a b -- ^ b-	-> Array a b -- ^ c+mv_mult :: (Ix i, Ix j, Num a) => Array (i,j) a -- ^ A+	-> Array j a -- ^ b+	-> Array i a -- ^ c -mv_mult a b = if ac /= br+mv_mult a b = if (ac0,ac1) /= bounds b 	      then error "mv_mult: dimensions inconsistent"-	      else array bnds [ (i, mult i) | i <- range bnds ]-    where mult i = sum [ a!(i,k) * b!(k) | k <- [1..ac] ]-	  ((_,_),(ar,ac)) = bounds a-	  (_,br) = bounds b-	  bnds = (1,ar)+	      else generate (ar0,ar1) $ \i ->+			sum [ a!(i,k) * bk | (k,bk) <- assocs b ]+    where ((ar0,ac0),(ar1,ac1)) = bounds a  -- | Transpose of a matrix -m_trans :: (Ix a, Integral a, Num b) => Array (a,a) b -- ^ A-	-> Array (a,a) b -- ^ A^T+m_trans :: (Ix i, Ix j, Num a) => Array (i,j) a -- ^ A+	-> Array (j,i) a -- ^ A^T -m_trans a = array bnds [ ((i,j), a!(j,i)) | (i,j) <- range bnds ]-    where (_,(m,n)) = bounds a-	  bnds = ((1,1),(n,m))+m_trans a = generate ((n0,m0),(n1,m1)) $ \(i,j) -> a!(j,i)+    where ((m0,n0),(m1,n1)) = bounds a  -- | Hermitian transpose (conjugate transpose) of a matrix -m_hermit :: (Ix a, Integral a, RealFloat b) => Array (a,a) (Complex b) -- ^ A-	 -> Array (a,a) (Complex b) -- ^ A^H+m_hermit :: (Ix i, Ix j, RealFloat a) => Array (i,j) (Complex a) -- ^ A+	 -> Array (j,i) (Complex a) -- ^ A^H -m_hermit a = array bnds [ ((i,j), conjugate (a!(j,i))) | (i,j) <- range bnds ]-    where (_,(m,n)) = bounds a-	  bnds = ((1,1),(n,m))+m_hermit a = generate ((n0,m0),(n1,m1)) $ \(i,j) -> conjugate (a!(j,i))+    where ((m0,n0),(m1,n1)) = bounds a+++columnBounds :: (Ix i, Ix j) => Array (i,j) a -> (i,i)+columnBounds a =+   let ((m0,_n0), (m1,_n1)) = bounds a+   in  (m0,m1)++rowBounds :: (Ix i, Ix j) => Array (i,j) a -> (j,j)+rowBounds a =+   let ((_m0,n0), (_m1,n1)) = bounds a+   in  (n0,n1)++getColumn :: (Ix i, Ix j) => j -> Array (i,j) e -> Array i e+getColumn j a = ixmap (columnBounds a) (\k -> (k,j)) a++getRow :: (Ix i, Ix j) => i -> Array (i,j) e -> Array j e+getRow k a = ixmap (rowBounds a) (\j -> (k,j)) a++toColumns :: (Ix i, Ix j) => Array (i,j) a -> [Array i a]+toColumns a = map (flip getColumn a) $ range $ rowBounds a++toRows :: (Ix i, Ix j) => Array (i,j) a -> [Array j a]+toRows a = map (flip getRow a) $ range $ columnBounds a+++{- |+We need the bounds of the row indices for empty input lists.+-}+fromColumns :: (Ix i) => (i,i) -> [Array i a] -> Array (i,Int) a+fromColumns bnds@(m0,m1) columns =+   if all ((bnds==) . bounds) columns+     then array ((m0,0), (m1, length columns - 1)) $ concat $+          zipWith+            (\k -> map (\(i,a) -> ((i,k),a)) . assocs)+            [0..] columns+     else error "Matrix.fromColumns: column bounds mismatch"++fromRows :: (Ix j) => (j,j) -> [Array j a] -> Array (Int,j) a+fromRows bnds@(n0,n1) rows =+   if all ((bnds==) . bounds) rows+     then array ((0,n0), (length rows - 1, n1)) $ concat $+          zipWith+            (\k -> map (\(i,a) -> ((k,i),a)) . assocs)+            [0..] rows+     else error "Matrix.fromRows: row bounds mismatch"++++outer :: (Ix i, Ix j, Num a) => Array i a -> Array j a -> Array (i,j) a+outer x y =+   let (m0,m1) = bounds x+       (n0,n1) = bounds y+   in  array ((m0,n0), (m1,n1)) $ do+         (i,xi) <- assocs x+         (j,yj) <- assocs y+         return ((i,j), xi*yj)++inner :: (Ix i, Num a) => Array i a -> Array i a -> a+inner x y =+   if bounds x == bounds y+     then sum $ zipWith (*) (elems x) (elems y)+     else error "inner: dimensions mismatch"

@@ -0,0 +1,160 @@+module Matrix.QR.Givens  (+   leastSquares,+   decompose, solve, det,+   Rotation, rotateVector,+   Upper, solveUpper, detUpper,+   ) where++import qualified Matrix.Sparse as Sparse+import DSP.Basic (toMaybe, (^!))++import Control.Monad (mfilter)++import qualified Data.Foldable as Fold+import qualified Data.List as List+import qualified Data.Map as Map+import Data.Map (Map)+import Data.Array+         (Array, Ix, array, bounds, elems, rangeSize, range, (!), (//), )+++data Rotation i a = Rotation (i,i) (a,a)+   deriving Show++data Upper i j a = Upper ((i,j), (i,j)) (Map i (Map j a))+   deriving Show++{- |+The decomposition routine is pretty simple.+It does not try to minimize fill-up by a clever ordering of rotations.+However, for banded matrices it will work as expected.+-}+decompose :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>+      Sparse.Matrix i j a -- ^ A+   -> ([Rotation i a], Upper i j a) -- ^ QR(A)+decompose a =+   (\(qs,rows) -> (concat qs, Upper (Sparse.bounds a) $ Map.fromList rows)) .+   unzip .+   List.unfoldr+      (\a0 ->+         let bnds@((m0,_), _) = Sparse.bounds a0+             (q,a1) = step a0+         in  toMaybe (not $ emptyRange bnds) $+               ((q, (m0, Sparse.getRow m0 a1)), submatrix a1))+    $ a++-- cf. QR.Householder+emptyRange :: (Ix i) => (i,i) -> Bool+emptyRange = null . range++-- | assumes that the first column is empty except the top-most element+submatrix ::+   (Ord i, Enum i, Ord j, Enum j) => Sparse.Matrix i j a -> Sparse.Matrix i j a+submatrix a =+   let ((m0,n0), mn1) = Sparse.bounds a+   in  Sparse.fromRows ((succ m0, succ n0), mn1) $+       Map.delete m0 $ Sparse.toRows a++step ::+   (Ord i, Ord j, RealFloat a) =>+   Sparse.Matrix i j a -> ([Rotation i a], Sparse.Matrix i j a)+step a =+   let bnds@((m0,n0), _) = Sparse.bounds a+       rows = Sparse.toRows a+       topRow = Map.findWithDefault Map.empty m0 rows+   in  (\((xi,xrem), finalRows) ->+         (Map.elems $ Map.mapMaybe fst finalRows,+          Sparse.fromRows bnds $ Map.insert m0 (Map.insert n0 xi xrem) $+            fmap snd finalRows)) $+       Map.mapAccumWithKey+         (\(xi,xrem) mi yrow ->+            let yrem = Map.delete n0 yrow+                jrot = Just . Rotation (m0,mi)+            in  case mfilter (0/=) $ Map.lookup n0 yrow of+                  Nothing -> ((xi,xrem),(Nothing,yrem))+                  Just yi ->+                     if xi==0+                        then ((yi,yrem), (jrot (0,-1), fmap negate xrem))+                        else+                           let rot = rotationFromXY (xi,yi)+                               (rx,ry) = rotateRows rot (xrem, yrem)+                           in  ((fst $ rotateSingle rot (xi,yi), rx),+                                (jrot rot, ry)))+         (Map.findWithDefault 0 n0 topRow, Map.delete n0 topRow)+         (Map.delete m0 rows)++-- | The argument must not be (0,0).+rotationFromXY :: (RealFloat a) => (a,a) -> (a,a)+rotationFromXY (x,y) =+   if abs x > abs y+     then let q = y/x; k = recipNorm q in (k,-q*k)+     else let q = x/y; k = recipNorm q in (q*k,-k)++recipNorm :: Floating a => a -> a+recipNorm q = recip $ sqrt (1+q^!2)++rotateSingle :: (Num a) => (a,a) -> (a,a) -> (a,a)+rotateSingle (c,s) (x,y) = (c*x-s*y, s*x+c*y)++rotateRows ::+   (Ord j, Num a) => (a,a) -> (Map j a, Map j a) -> (Map j a, Map j a)+rotateRows (c,s) (xs,ys) =+   let rs =+         Map.intersectionWith (curry $ rotateSingle (c,s)) xs ys+         `Map.union`+         fmap (\x -> ( c*x, s*x)) (Map.difference xs ys)+         `Map.union`+         fmap (\y -> (-s*y, c*y)) (Map.difference ys xs)++   in  (fmap fst rs, fmap snd rs)++rotateVector :: (Ix i, Num a) => Rotation i a -> Array i a -> Array i a+rotateVector (Rotation (i0,i1) cs) x =+   let (y0,y1) = rotateSingle cs (x!i0,x!i1)+   in  x // [(i0,y0),(i1,y1)]+++{- |+Assumes that 'Upper' matrix is at least as high as wide+and that it has full rank.+-}+solveUpper ::+   (Ix i, Ix j, Fractional a) => Upper i j a -> Array i a -> Array j a+solveUpper (Upper ((m0,n0), (m1,n1)) rows0) b =+   if bounds b == (m0,m1)+     then+         array (n0,n1) $ Map.toList $+         foldr+            (\(row,bi) xs ->+               let ((j,a),as) = Map.deleteFindMin row+               in  Map.insert j+                     ((bi - Fold.sum (Map.intersectionWith (*) as xs)) / a) xs)+            Map.empty+            (zip (Map.elems rows0) (elems b))+     else error "solveUpper: vertical bounds mismatch"++solve ::+   (Ix i, Ix j, Fractional a) =>+   ([Rotation i a], Upper i j a) -> Array i a -> Array j a+solve (qs, u) b = solveUpper u $ foldl (flip rotateVector) b qs++{- |+Solve a sparse overconstrained linear problem, i.e. minimize @||Ax-b||@.+@A@ must have dimensions @m x n@ with @m>=n@ and it must have full-rank.+None of these conditions is checked.+-}+leastSquares ::+   (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>+   Sparse.Matrix i j a -> Array i a -> Array j a+leastSquares = solve . decompose+++detUpper ::+   (Ix i, Ix j, Fractional a) => Upper i j a -> a+detUpper (Upper ((_m0,n0), (_m1,n1)) rows) =+   if rangeSize (n0,n1) == Map.size rows+     then product $ map (snd . Map.findMin) $ Map.elems rows+     else 0++det :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Sparse.Matrix i j a -> a+det = detUpper . snd . decompose

@@ -0,0 +1,152 @@+module Matrix.QR.Householder (+   leastSquares,+   decompose, solve, det,+   Reflection, reflectMatrix, reflectVector,+   Upper, matrixFromUpper, solveUpper, detUpper,+   ) where++import Matrix.Matrix (mv_mult, m_trans, getRow, getColumn, inner, outer)+import Matrix.Vector (sub, scale, norm)+import DSP.Basic (toMaybe)++import qualified Data.List as List+import Data.Array+         (Array, Ix, bounds, elems, range, rangeSize,+          accum, accumArray, assocs, ixmap, listArray, (!), (//), )+++decompose :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>+      Array (i,j) a -- ^ A+   -> ([Reflection i a], Upper i j a) -- ^ QR(A)+decompose a =+   (\(qs,rows) -> (qs, Upper (bounds a) rows)) .+   unzip .+   List.unfoldr+      (\a0 ->+         let bnds@((m0,_), _) = bounds a0+         in  toMaybe (not $ emptyRange bnds) $+               let (q,a1) = step a0+               in  ((q, getRow m0 a1), submatrix a1))+    $ a++emptyRange :: (Ix i) => (i,i) -> Bool+emptyRange = null . range++step ::+   (Ix i, Ix j, RealFloat a) =>+   Array (i,j) a -> (Reflection i a, Array (i,j) a)+step a =+   let (m0,n0) = fst $ bounds a+       z = getColumn n0 a+       sign x = if x<0 then -1 else 1+       q = reflection $ accum (+) z [(m0, sign(z!m0) * norm z)]+   in  (q, reflectMatrix q a)++{-+Submatrices with only Ix constrained indices would not work,+because we cannot reduce a two-dimensional array by only one element.+-}+submatrix :: (Ix i, Enum i, Ix j, Enum j) => Array (i,j) e -> Array (i,j) e+submatrix a =+   let ((m0,n0), (m1,n1)) = bounds a+   in  ixmap ((succ m0, succ n0), (m1,n1)) id a+++data Upper i j a = Upper ((i,j), (i,j)) [Array j a]++matrixFromUpper :: (Ix i, Ix j, Num a) => Upper i j a -> Array (i,j) a+matrixFromUpper (Upper bnds@((m0,_n0), (m1,_n1)) rows) =+   accumArray (const id) 0 bnds $ concat $+   zipWith (\k -> map (\(j,a) -> ((k,j),a)) . assocs) (range (m0,m1)) rows+++newtype Reflection i a = Reflection (Array i a)++reflection :: (Ix i, Floating a) => Array i a -> Reflection i a+reflection v =+   let normv = norm v+   in  Reflection $ fmap (/ ((1-signum normv) + normv)) v++reflectMatrixFull ::+   (Ix i, Ix j, Num a) => Reflection i a -> Array (i,j) a -> Array (i,j) a+reflectMatrixFull (Reflection v) a =+   sub a $ scale 2 $ outer v $ mv_mult (m_trans a) v++reflectMatrix ::+   (Ix i, Ix j, Num a) => Reflection i a -> Array (i,j) a -> Array (i,j) a+reflectMatrix q@(Reflection v) a =+   let (k0,k1) = bounds v+       ((m0,n0), (m1,n1)) = bounds a+       bnds = ((k0,n0),(k1,n1))+   in  case (compare k0 m0, compare k1 m1) of+         (EQ,EQ) -> reflectMatrixFull q a+         (LT,_) -> error "reflectMatrix: lower reflection dimension too small"+         (_,GT) -> error "reflectMatrix: upper reflection dimension too big"+         _ -> replaceSubArray a $ reflectMatrixFull q $ subArray bnds a+++reflectVectorFull :: (Ix i, Num a) => Reflection i a -> Array i a -> Array i a+reflectVectorFull (Reflection v) a = sub a $ scale (2 * inner v a) v++reflectVector :: (Ix i, Num a) => Reflection i a -> Array i a -> Array i a+reflectVector q@(Reflection v) a =+   let bnds@(k0,k1) = bounds v+       (m0,m1) = bounds a+   in  case (compare k0 m0, compare k1 m1) of+         (EQ,EQ) -> reflectVectorFull q a+         (LT,_) -> error "reflectVector: lower reflection dimension too small"+         (_,GT) -> error "reflectVector: upper reflection dimension too big"+         _ -> replaceSubArray a $ reflectVectorFull q $ subArray bnds a+++subArray :: (Ix i) => (i,i) -> Array i a -> Array i a+subArray bnds = ixmap bnds id++replaceSubArray :: (Ix i) => Array i a -> Array i a -> Array i a+replaceSubArray x y = x // assocs y+++{- |+Assumes that 'Upper' matrix is at least as high as wide+and that it has full rank.+-}+solveUpper ::+   (Ix i, Ix j, Fractional a) => Upper i j a -> Array i a -> Array j a+solveUpper (Upper ((m0,n0), (m1,n1)) rs0) b =+   if bounds b == (m0,m1)+     then+         listArray (n0,n1) $+         foldr+            (\(r,bi) xs ->+               let (a:as) = elems r+               in  (bi - sum (zipWith (*) as xs)) / a : xs)+            []+            (zip rs0 (elems b))+     else error "solveUpper: vertical bounds mismatch"++solve ::+   (Ix i, Ix j, Fractional a) =>+   ([Reflection i a], Upper i j a) -> Array i a -> Array j a+solve (qs, u) b = solveUpper u $ foldl (flip reflectVector) b qs++{- |+Solve an overconstrained linear problem, i.e. minimize @||Ax-b||@.+@A@ must have dimensions @m x n@ with @m>=n@ and it must have full-rank.+None of these conditions is checked.+-}+leastSquares ::+   (Ix i, Enum i, Ix j, Enum j, RealFloat a) =>+   Array (i,j) a -> Array i a -> Array j a+leastSquares = solve . decompose+++detUpper :: (Ix i, Ix j, Fractional a) => Upper i j a -> a+detUpper (Upper ((_m0,n0), (_m1,n1)) rs) =+   if rangeSize (n0,n1) == length rs+     then product $ map (head . elems) rs+     else 0++det :: (Ix i, Enum i, Ix j, Enum j, RealFloat a) => Array (i,j) a -> a+det a =+   let (qs,u) = decompose a+   in  (if even (length qs) then 1 else -1) * detUpper u

@@ -0,0 +1,84 @@+module Matrix.Sparse (+   Matrix,+   bounds,+   fromMap,+   fromRows,+   fromColumns,+   fromDense,+   toRows,+   toColumns,+   toDense,+   getRow,+   getColumn,+   mulVector,+   ) where++import qualified Matrix.Vector as Vector+import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import qualified Data.Array as Array+import Data.Map (Map)+import Data.Array (Array, Ix, accumArray, (!))+++data Matrix i j a = Matrix ((i,j), (i,j)) (Map i (Map j a))+   deriving Show++instance Functor (Matrix i j) where+   fmap f (Matrix bnds m) = Matrix bnds $ fmap (fmap f) m+++bounds :: Matrix i j a -> ((i,j), (i,j))+bounds (Matrix bnds _) = bnds++fromMap :: (Ord i, Ord j) => ((i,j), (i,j)) -> Map (i,j) a -> Matrix i j a+fromMap bnds =+   Matrix bnds . Map.fromListWith Map.union .+   map (\((i,j),a) -> (i, Map.singleton j a)) . Map.toList++fromRows ::+   (Ord i, Ord j) => ((i,j), (i,j)) -> Map i (Map j a) -> Matrix i j a+fromRows = Matrix++fromColumns ::+   (Ord i, Ord j) => ((i,j), (i,j)) -> Map j (Map i a) -> Matrix i j a+fromColumns bnds = Matrix bnds . flipMap++fromDense :: (Ix i, Ix j) => Array (i,j) a -> Matrix i j a+fromDense a = fromMap (Array.bounds a) $ Map.fromList $ Array.assocs a+++toRows :: (Ord i, Ord j) => Matrix i j a -> Map i (Map j a)+toRows (Matrix _bnds rows) = rows++toColumns :: (Ord i, Ord j) => Matrix i j a -> Map j (Map i a)+toColumns (Matrix _bnds rows) = flipMap rows++toDense :: (Ix i, Ix j, Num a) => Matrix i j a -> Array (i,j) a+toDense (Matrix bnds a) =+   accumArray (const id) 0 bnds $ Fold.fold $+   Map.mapWithKey (\i -> map (\(j,e) -> ((i,j),e)) .  Map.toList) a+++-- cf. comfort-graph:Graph.Comfort.Map.flip+flipMap :: (Ord i, Ord j) => Map i (Map j a) -> Map j (Map i a)+flipMap =+   Map.unionsWith (Map.unionWith (error $ "Map.flip: duplicate key")) .+   Map.elems . Map.mapWithKey (fmap . Map.singleton)+++getRow :: (Ord i, Ord j) => i -> Matrix i j a -> Map j a+getRow i (Matrix _ rows) = Map.findWithDefault Map.empty i rows++getColumn :: (Ord i, Ord j) => j -> Matrix i j a -> Map i a+getColumn j (Matrix _ rows) = Map.mapMaybe (Map.lookup j) rows+++mulVector :: (Ix i, Ix j, Num a) => Matrix i j a -> Array j a -> Array i a+mulVector a@(Matrix ((m0,n0), (m1,n1)) _) v =+   if (n0,n1) == Array.bounds v+     then Vector.generate (m0,m1) $ flip mulRowVector v . flip getRow a+     else error "Sparse.mulVector: dimensions mismatch"++mulRowVector :: (Ix j, Num a) => Map j a -> Array j a -> a+mulRowVector row v = Fold.sum $ Map.mapWithKey (\j x -> x * v!j) row

@@ -0,0 +1,35 @@+module Matrix.Vector where++import DSP.Basic ((^!))+import Data.Array+         (Array, Ix, bounds, elems, range, array, assocs, listArray, (!), )+++generate :: (Ix i) => (i,i) -> (i -> a) -> Array i a+generate bnds f = array bnds $ map (\i -> (i, f i)) $ range bnds++fromList :: [a] -> Array Int a+fromList xs = listArray (0, length xs - 1) xs++toList :: Array Int a -> [a]+toList = elems+++norm :: (Ix i, Floating a) => Array i a -> a+norm = sqrt . sum . elems . fmap (^!2)++scale :: (Ix i, Num a) => a -> Array i a -> Array i a+scale x = fmap (x*)+++lift2 :: (Ix i) => (a -> b -> c) -> Array i a -> Array i b -> Array i c+lift2 f x y =+   if bounds x == bounds y+     then array (bounds x) [ (k, f xk (y!k)) | (k, xk) <- assocs x ]+     else error "Vector.lift2: matrix dimensions mismatch"++add :: (Ix i, Num a) => Array i a -> Array i a -> Array i a+add = lift2 (+)++sub :: (Ix i, Num a) => Array i a -> Array i a -> Array i a+sub = lift2 (-)

@@ -1,5 +1,5 @@ Name:             dsp-Version:          0.2.3.1+Version:          0.2.4 License:          GPL License-File:     COPYING Copyright:        Matt Donadio, 2003@@ -11,14 +11,13 @@ Synopsis:         Haskell Digital Signal Processing Description:      Digital Signal Processing, Fourier Transform, Linear Algebra, Interpolation Category:         Sound, Math--Tested-With:      GHC==6.4.1, GHC==6.8.2-Tested-With:      GHC==7.4.2, GHC==7.6.3+Tested-With:      GHC==7.4.2, GHC==7.6.3, GHC==7.8.4+Tested-With:      GHC==8.0.2, GHC==8.2.1 Cabal-Version:    >=1.8 Build-Type:       Simple  Source-Repository this-  Tag:         0.2.3.1+  Tag:         0.2.4   Type:        darcs   Location:    http://code.haskell.org/~thielema/dsp/ @@ -26,22 +25,16 @@   Type:        darcs   Location:    http://code.haskell.org/~thielema/dsp/ -Flag splitBase-  Description: Choose the new smaller, split-up base package.- Flag buildExamples   Description: Build demo executables-  Default: True+  Default: False  Library-  If flag(splitBase)-    Build-Depends:-      array >=0.1 && <0.6,-      random >=1.0 && <1.2,-      base >= 2 && <5-  Else-    Build-Depends:-      base >=1.0 && <2+  Build-Depends:+    containers >=0.3 && <0.6,+    array >=0.1 && <0.6,+    random >=1.0 && <1.2,+    base >= 2 && <5    GHC-Options:      -Wall   Exposed-modules:@@ -75,6 +68,7 @@     DSP.Filter.IIR.Matchedz     DSP.Filter.IIR.Prony     DSP.Filter.IIR.Transform+    DSP.Filter.IIR.Cookbook     DSP.Flowgraph     DSP.Multirate.CIC     DSP.Multirate.Halfband@@ -85,8 +79,12 @@     DSP.Window     Matrix.Cholesky     Matrix.LU+    Matrix.QR.Householder+    Matrix.QR.Givens     Matrix.Levinson     Matrix.Matrix+    Matrix.Vector+    Matrix.Sparse     Matrix.Simplex     Numeric.Approximation.Chebyshev     Numeric.Random.Distribution.Binomial@@ -123,9 +121,21 @@     Polynomial.Chebyshev     Polynomial.Maclaurin     Polynomial.Roots-    DSP.Filter.IIR.Cookbook   Other-Modules:     Numeric.Transform.Fourier.Eigensystem++Test-Suite dsp-test+  Type: exitcode-stdio-1.0+  Main-Is: Main.hs+  Other-Modules: Test.Matrix.QR+  Hs-Source-Dirs: test+  GHC-Options: -Wall+  Build-Depends:+    QuickCheck >=2.5 && <3,+    dsp,+    containers,+    array,+    base  Executable dsp-demo-article   Main-Is: Article.hs

@@ -0,0 +1,15 @@+module Main where++import qualified Test.Matrix.QR as QR+++prefix :: String -> [(String, IO ())] -> [(String, IO ())]+prefix msg =+   map (\(str,test) -> (msg ++ "." ++ str, test))++main :: IO ()+main =+   mapM_ (\(msg,io) -> putStr (msg++": ") >> io) $+   concat $+      prefix "Matrix.QR" QR.tests :+      []

@@ -0,0 +1,164 @@+module Test.Matrix.QR where++import qualified Matrix.QR.Householder as Householder+import qualified Matrix.QR.Givens as Givens+import qualified Matrix.LU as LU+import qualified Matrix.Vector as Vector+import qualified Matrix.Sparse as Sparse+import Matrix.Matrix (m_trans, mm_mult, mv_mult)+import DSP.Basic ((^!))++import Control.Applicative (liftA2, (<$>))++import qualified Data.Foldable as Fold+import qualified Data.Map as Map+import Data.Array++import qualified Test.QuickCheck as QC+++doubleArray :: (Ix i) => Array i Int -> Array i Double+doubleArray = fmap fromIntegral++gramian :: Array (Int,Int) Double -> Double+gramian m = LU.det (m_trans m `mm_mult` m)++fullRank :: Array (Int,Int) Int -> Bool+fullRank m = round (gramian $ doubleArray m) /= (0::Integer)++arbitraryIntArray :: (Ix i) => (i,i) -> QC.Gen (Array i Int)+arbitraryIntArray bnds =+   fmap (listArray bnds) $ QC.vectorOf (rangeSize bnds) $ QC.choose (-10,10)++genMatrix :: QC.Gen (Array (Int,Int) Int)+genMatrix = do+   m <- QC.choose (0,5)+   n <- QC.choose (0,m)+   arbitraryIntArray ((1,1),(m,n))++genForward :: QC.Gen (Array (Int,Int) Int, Array Int Int)+genForward = do+   a <- genMatrix `QC.suchThat` fullRank+   let ((_m0,n0), (_m1,n1)) = bounds a+   x <- arbitraryIntArray (n0,n1)+   return (a,x)++genInverse :: QC.Gen (Array (Int,Int) Int, Array Int Int)+genInverse = do+   a <- genMatrix `QC.suchThat` fullRank+   let ((m0,_n0), (m1,_n1)) = bounds a+   b <- arbitraryIntArray (m0,m1)+   return (a,b)+++arbitraryIntSparse ::+   (Ix i, Ix j) => ((i, j), (i, j)) -> QC.Gen (Sparse.Matrix i j Int)+arbitraryIntSparse bnds =+   fmap+      (Sparse.fromMap bnds . fmap fst .+       Map.filter snd . Map.fromList . zip (range bnds)) $+      QC.vectorOf (rangeSize bnds) $+      liftA2 (,) (QC.choose (-10,10)) QC.arbitrary++genSparse :: QC.Gen (Sparse.Matrix Int Int Int)+genSparse = do+   m <- QC.choose (0,5)+   n <- QC.choose (0,m)+   arbitraryIntSparse ((1,1),(m,n))++genSparseInverse :: QC.Gen (Sparse.Matrix Int Int Int, Array Int Int)+genSparseInverse = do+   a <- genSparse `QC.suchThat` (fullRank . Sparse.toDense)+   let ((m0,_n0), (m1,_n1)) = Sparse.bounds a+   b <- arbitraryIntArray (m0,m1)+   return (a,b)++genSquare :: QC.Gen (Array (Int,Int) Int)+genSquare = do+   m <- QC.choose (0,5)+   arbitraryIntArray ((1,1),(m,m))++genSparseSquare :: QC.Gen (Sparse.Matrix Int Int Int)+genSparseSquare = do+   m <- QC.choose (0,5)+   arbitraryIntSparse ((1,1),(m,m))++++approx :: (Fractional a, Ord a) => a -> a -> Bool+approx x y  =  abs (x-y) <= 1e-5 * max 1 (abs x + abs y)++approxAbsVector :: (Fractional a, Ord a, Ix i) => Array i a -> Array i a -> Bool+approxAbsVector x y = (Fold.foldl max 0 $ fmap abs $ Vector.sub x y) < 1e-5++++solveHouseholder :: QC.Property+solveHouseholder =+   QC.forAll genForward $ \(a,x) ->+      let b = mv_mult a x+      in  x ==+          fmap round (Householder.leastSquares (doubleArray a) (doubleArray b))++solveGivens :: QC.Property+solveGivens =+   QC.forAll genForward $ \(a,x) ->+      let b = mv_mult a x+      in  x ==+          fmap round+            (Givens.leastSquares+               (Sparse.fromDense $ doubleArray a) (doubleArray b))+++leastSquares :: QC.Property+leastSquares =+   QC.forAll genInverse $ \(a,b) ->+      Householder.leastSquares (doubleArray a) (doubleArray b)+      `approxAbsVector`+      Givens.leastSquares (Sparse.fromDense $ doubleArray a) (doubleArray b)++leastSquaresSparse :: QC.Property+leastSquaresSparse =+   QC.forAll genSparseInverse $ \(a,b) ->+      Householder.leastSquares (doubleArray $ Sparse.toDense a) (doubleArray b)+      `approxAbsVector`+      Givens.leastSquares (fmap fromIntegral a) (doubleArray b)++++gramianHouseholder :: QC.Property+gramianHouseholder =+   QC.forAll (fmap doubleArray genMatrix) $ \a ->+      gramian a `approx` (Householder.det a ^! 2)++gramianGivens :: QC.Property+gramianGivens =+   QC.forAll (fmap fromIntegral <$> genSparse) $ \a ->+      gramian (Sparse.toDense a)  `approx`  (Givens.det a ^! 2)++detHouseholder :: QC.Property+detHouseholder =+   QC.forAll (fmap doubleArray genSquare) $ \a ->+      LU.det a `approx` Householder.det a++detGivens :: QC.Property+detGivens =+   QC.forAll (fmap fromIntegral <$> genSparseSquare) $ \a ->+      LU.det (Sparse.toDense a)  `approx`  Givens.det a+++longCheck :: QC.Property -> IO ()+longCheck =+   QC.quickCheckWith (QC.stdArgs {QC.maxSuccess=10000})++tests :: [(String, IO ())]+tests =+   ("solveHouseholder", longCheck solveHouseholder) :+   ("solveGivens", longCheck solveGivens) :+   ("leastSquares", longCheck leastSquares) :+   ("leastSquaresSparse", longCheck leastSquaresSparse) :+   ("gramianHouseholder", longCheck gramianHouseholder) :+   ("gramianGivens", longCheck gramianGivens) :+   ("detHouseholder", longCheck detHouseholder) :+   ("detGivens", longCheck detGivens) :+   []