X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=src%2FLinear%2FMatrix.hs;h=39576dc41343c2bbad77b745a33bc720d44f666f;hb=303c5e7bba583f08e59bc6c848be8e75c1155a3b;hp=63c0348cd9e6e0ce676ce5c17c1629890836a6e6;hpb=d52e10c90c0b8263af2e6a0152cebf0ad3c70e62;p=numerical-analysis.git diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index 63c0348..39576dc 100644 --- a/src/Linear/Matrix.hs +++ b/src/Linear/Matrix.hs @@ -1,78 +1,146 @@ -{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ExistentialQuantification #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE RebindableSyntax #-} module Linear.Matrix where +import Data.List (intercalate) + import Data.Vector.Fixed ( Dim, Vector ) import qualified Data.Vector.Fixed as V ( + Fun(..), + N1, + and, + eq, + foldl, fromList, length, map, - toList + maximum, + replicate, + toList, + zipWith ) -import Data.Vector.Fixed.Internal (arity) - +import Data.Vector.Fixed.Internal (Arity, arity, S, Dim) import Linear.Vector +import Normed + +import NumericPrelude hiding (abs) +import qualified Algebra.Algebraic as Algebraic +import qualified Algebra.Absolute as Absolute +import qualified Algebra.Additive as Additive +import qualified Algebra.Ring as Ring +import Algebra.Absolute (abs) +import qualified Algebra.Field as Field +import qualified Algebra.RealField as RealField +import qualified Algebra.RealRing as RealRing +import qualified Algebra.ToRational as ToRational +import qualified Algebra.Transcendental as Transcendental +import qualified Prelude as P + +data Mat v w a = (Vector v (w a), Vector w a) => Mat (v (w a)) +type Mat1 a = Mat D1 D1 a +type Mat2 a = Mat D2 D2 a +type Mat3 a = Mat D3 D3 a +type Mat4 a = Mat D4 D4 a + +-- We can't just declare that all instances of Vector are instances of +-- Eq unfortunately. We wind up with an overlapping instance for +-- w (w a). +instance (Eq a, Vector v Bool, Vector w Bool) => Eq (Mat v w a) where + -- | Compare a row at a time. + -- + -- Examples: + -- + -- >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int + -- >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int + -- >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int + -- >>> m1 == m2 + -- True + -- >>> m1 == m3 + -- False + -- + (Mat rows1) == (Mat rows2) = + V.and $ V.zipWith comp rows1 rows2 + where + -- Compare a row, one column at a time. + comp row1 row2 = V.and (V.zipWith (==) row1 row2) + + +instance (Show a, Vector v String, Vector w String) => Show (Mat v w a) where + -- | Display matrices and vectors as ordinary tuples. This is poor + -- practice, but these results are primarily displayed + -- interactively and convenience trumps correctness (said the guy + -- who insists his vector lengths be statically checked at + -- compile-time). + -- + -- Examples: + -- + -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int + -- >>> show m + -- ((1,2),(3,4)) + -- + show (Mat rows) = + "(" ++ (intercalate "," (V.toList row_strings)) ++ ")" + where + row_strings = V.map show_vector rows + show_vector v1 = + "(" ++ (intercalate "," element_strings) ++ ")" + where + v1l = V.toList v1 + element_strings = P.map show v1l + -type Mat v w a = Vn v (Vn w a) -type Mat2 a = Mat Vec2D Vec2D a -type Mat3 a = Mat Vec3D Vec3D a -type Mat4 a = Mat Vec4D Vec4D a -- | Convert a matrix to a nested list. -toList :: (Vector v (Vn w a), Vector w a) => Mat v w a -> [[a]] -toList m = map V.toList (V.toList m) +toList :: Mat v w a -> [[a]] +toList (Mat rows) = map V.toList (V.toList rows) -- | Create a matrix from a nested list. -fromList :: (Vector v (Vn w a), Vector w a) => [[a]] -> Mat v w a -fromList vs = V.fromList $ map V.fromList vs +fromList :: (Vector v (w a), Vector w a, Vector v a) => [[a]] -> Mat v w a +fromList vs = Mat (V.fromList $ map V.fromList vs) -- | Unsafe indexing. -(!!!) :: (Vector v (Vn w a), Vector w a) => Mat v w a -> (Int, Int) -> a +(!!!) :: (Vector w a) => Mat v w a -> (Int, Int) -> a (!!!) m (i, j) = (row m i) ! j -- | Safe indexing. -(!!?) :: (Vector v (Vn w a), Vector w a) => Mat v w a - -> (Int, Int) - -> Maybe a -(!!?) m (i, j) +(!!?) :: Mat v w a -> (Int, Int) -> Maybe a +(!!?) m@(Mat rows) (i, j) | i < 0 || j < 0 = Nothing - | i > V.length m = Nothing + | i > V.length rows = Nothing | otherwise = if j > V.length (row m j) then Nothing else Just $ (row m j) ! j -- | The number of rows in the matrix. -nrows :: forall v w a. (Vector v (Vn w a), Vector w a) => Mat v w a -> Int -nrows = V.length +nrows :: Mat v w a -> Int +nrows (Mat rows) = V.length rows -- | The number of columns in the first row of the -- matrix. Implementation stolen from Data.Vector.Fixed.length. -ncols :: forall v w a. (Vector v (Vn w a), Vector w a) => Mat v w a -> Int -ncols _ = arity (undefined :: Dim w) +ncols :: forall v w a. (Vector w a) => Mat v w a -> Int +ncols _ = (arity (undefined :: Dim w)) -- | Return the @i@th row of @m@. Unsafe. -row :: (Vector v (Vn w a), Vector w a) => Mat v w a - -> Int - -> Vn w a -row m i = m ! i +row :: Mat v w a -> Int -> w a +row (Mat rows) i = rows ! i -- | Return the @j@th column of @m@. Unsafe. -column :: (Vector v a, Vector v (Vn w a), Vector w a) => Mat v w a - -> Int - -> Vn v a -column m j = - V.map (element j) m +column :: (Vector v a) => Mat v w a -> Int -> v a +column (Mat rows) j = + V.map (element j) rows where element = flip (!) @@ -89,16 +157,16 @@ column m j = -- >>> transpose m -- ((1,3),(2,4)) -- -transpose :: (Vector v (Vn w a), - Vector w (Vn v a), +transpose :: (Vector w (v a), Vector v a, Vector w a) => Mat v w a -> Mat w v a -transpose m = V.fromList column_list +transpose m = Mat $ V.fromList column_list where column_list = [ column m i | i <- [0..(ncols m)-1] ] + -- | Is @m@ symmetric? -- -- Examples: @@ -111,7 +179,7 @@ transpose m = V.fromList column_list -- >>> symmetric m2 -- False -- -symmetric :: (Vector v (Vn w a), +symmetric :: (Vector v (w a), Vector w a, v ~ w, Vector w Bool, @@ -136,11 +204,11 @@ symmetric m = -- ((0,1,2),(1,2,3),(2,3,4)) -- construct :: forall v w a. - (Vector v (Vn w a), + (Vector v (w a), Vector w a) => (Int -> Int -> a) -> Mat v w a -construct lambda = rows +construct lambda = Mat rows where -- The arity trick is used in Data.Vector.Fixed.length. imax = (arity (undefined :: Dim v)) - 1 @@ -161,19 +229,21 @@ construct lambda = rows -- ((20.000000000000004,-1.0),(-1.0,20.0)) -- cholesky :: forall a v w. - (RealFloat a, - Vector v (Vn w a), - Vector w a) + (Algebraic.C a, + Vector v (w a), + Vector w a, + Vector v a) => (Mat v w a) -> (Mat v w a) cholesky m = construct r where r :: Int -> Int -> a - r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)**2 | k <- [0..i-1]]) + r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]]) | i < j = (((m !!! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i) | otherwise = 0 + -- | Matrix multiplication. Our 'Num' instance doesn't define one, and -- we need additional restrictions on the result type anyway. -- @@ -184,12 +254,11 @@ cholesky m = construct r -- >>> m1 `mult` m2 -- ((22,28),(49,64)) -- -mult :: (Num a, - Vector v (Vn w a), +mult :: (Ring.C a, + Vector v a, Vector w a, - Vector w (Vn z a), Vector z a, - Vector v (Vn z a)) + Vector v (z a)) => Mat v w a -> Mat w z a -> Mat v z a @@ -197,3 +266,104 @@ mult m1 m2 = construct lambda where lambda i j = sum [(m1 !!! (i,k)) * (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ] + + + +instance (Ring.C a, + Vector v (w a), + Vector w a) + => Additive.C (Mat v w a) where + + (Mat rows1) + (Mat rows2) = + Mat $ V.zipWith (V.zipWith (+)) rows1 rows2 + + (Mat rows1) - (Mat rows2) = + Mat $ V.zipWith (V.zipWith (-)) rows1 rows2 + + zero = Mat (V.replicate $ V.replicate (fromInteger 0)) + + +instance (Ring.C a, + Vector v (w a), + Vector w a, + v ~ w) + => Ring.C (Mat v w a) where + one = Mat (V.replicate $ V.replicate (fromInteger 1)) + m1 * m2 = m1 `mult` m2 + + +instance (Algebraic.C a, + ToRational.C a, + Vector v (w a), + Vector w a, + Vector v a, + Vector v [a]) + => Normed (Mat v w a) where + -- Treat the matrix as a big vector. + norm_p p (Mat rows) = + sqrt $ sum [(fromRational' $ toRational x)^2 | x <- xs] + where + xs = concat $ V.toList $ V.map V.toList rows + + norm_infty m@(Mat rows) + | nrows m == 0 || ncols m == 0 = 0 + | otherwise = + fromRational' $ toRational $ + P.maximum $ V.toList $ V.map (P.maximum . V.toList) rows + + + + + +-- Vector helpers. We want it to be easy to create low-dimension +-- column vectors. +type Vec a b = Mat a D1 b + +vec2d :: (a,a) -> Mat D2 D1 a +vec2d (x,y) = Mat (D2 (D1 x) (D1 y)) + +vec3d :: (a,a,a) -> Mat D3 D1 a +vec3d (x,y,z) = Mat (D3 (D1 x) (D1 y) (D1 z)) + +vec4d :: (a,a,a,a) -> Mat D4 D1 a +vec4d (w,x,y,z) = Mat (D4 (D1 w) (D1 x) (D1 y) (D1 z)) + +dot :: (RealRing.C a, + Dim w ~ V.N1, + Vector v a, + Vector w a, + Vector w (v a), + Vector w (w a)) + => Mat v w a + -> Mat v w a + -> a +v1 `dot` v2 = ((transpose v1) `mult` v2) !!! (0, 0) + + +-- | The angle between @v1@ and @v2@ in Euclidean space. +-- +-- Examples: +-- +-- >>> let v1 = make2d (1.0, 0.0) +-- >>> let v2 = make2d (0.0, 1.0) +-- >>> angle v1 v2 == pi/2.0 +-- True +-- +angle :: (Transcendental.C a, + RealRing.C a, + Dim w ~ V.N1, + Vector w (w a), + Vector v [a], + Vector v a, + Vector w a, + Vector v (w a), + Vector w (v a), + ToRational.C a) + => Mat v w a + -> Mat v w a + -> a +angle v1 v2 = + acos theta + where + theta = (recip norms) * (v1 `dot` v2) + norms = (norm v1) * (norm v2)