src/Linear/Matrix.hs

   1 {-# LANGUAGE ExistentialQuantification #-}
   2 {-# LANGUAGE FlexibleContexts #-}
   3 {-# LANGUAGE FlexibleInstances #-}
   4 {-# LANGUAGE MultiParamTypeClasses #-}
   5 {-# LANGUAGE ScopedTypeVariables #-}
   6 {-# LANGUAGE TypeFamilies #-}
   7 {-# LANGUAGE RebindableSyntax #-}
   8
   9 module Linear.Matrix
  10 where
  11
  12 import Data.List (intercalate)
  13
  14 import Data.Vector.Fixed (
  15   Dim,
  16   Vector
  17   )
  18 import qualified Data.Vector.Fixed as V (
  19   Fun(..),
  20   N1,
  21   and,
  22   eq,
  23   foldl,
  24   fromList,
  25   length,
  26   map,
  27   maximum,
  28   replicate,
  29   toList,
  30   zipWith
  31   )
  32 import Data.Vector.Fixed.Internal (Arity, arity, S, Dim)
  33 import Linear.Vector
  34 import Normed
  35
  36 import NumericPrelude hiding (abs)
  37 import qualified Algebra.Algebraic as Algebraic
  38 import qualified Algebra.Absolute as Absolute
  39 import qualified Algebra.Additive as Additive
  40 import qualified Algebra.Ring as Ring
  41 import Algebra.Absolute (abs)
  42 import qualified Algebra.Field as Field
  43 import qualified Algebra.RealField as RealField
  44 import qualified Algebra.RealRing as RealRing
  45 import qualified Algebra.ToRational as ToRational
  46 import qualified Algebra.Transcendental as Transcendental
  47 import qualified Prelude as P
  48
  49 data Mat v w a = (Vector v (w a), Vector w a) => Mat (v (w a))
  50 type Mat1 a = Mat D1 D1 a
  51 type Mat2 a = Mat D2 D2 a
  52 type Mat3 a = Mat D3 D3 a
  53 type Mat4 a = Mat D4 D4 a
  54
  55 -- We can't just declare that all instances of Vector are instances of
  56 -- Eq unfortunately. We wind up with an overlapping instance for
  57 -- w (w a).
  58 instance (Eq a, Vector v Bool, Vector w Bool) => Eq (Mat v w a) where
  59   -- | Compare a row at a time.
  60   --
  61   --   Examples:
  62   --
  63   --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
  64   --   >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
  65   --   >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
  66   --   >>> m1 == m2
  67   --   True
  68   --   >>> m1 == m3
  69   --   False
  70   --
  71   (Mat rows1) == (Mat rows2) =
  72     V.and $ V.zipWith comp rows1 rows2
  73     where
  74       -- Compare a row, one column at a time.
  75       comp row1 row2 = V.and (V.zipWith (==) row1 row2)
  76
  77
  78 instance (Show a, Vector v String, Vector w String) => Show (Mat v w a) where
  79   -- | Display matrices and vectors as ordinary tuples. This is poor
  80   --   practice, but these results are primarily displayed
  81   --   interactively and convenience trumps correctness (said the guy
  82   --   who insists his vector lengths be statically checked at
  83   --   compile-time).
  84   --
  85   --   Examples:
  86   --
  87   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
  88   --   >>> show m
  89   --   ((1,2),(3,4))
  90   --
  91   show (Mat rows) =
  92     "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
  93     where
  94       row_strings = V.map show_vector rows
  95       show_vector v1 =
  96         "(" ++ (intercalate "," element_strings) ++ ")"
  97         where
  98           v1l = V.toList v1
  99           element_strings = P.map show v1l
 100
 101
 102
 103 -- | Convert a matrix to a nested list.
 104 toList :: Mat v w a -> [[a]]
 105 toList (Mat rows) = map V.toList (V.toList rows)
 106
 107 -- | Create a matrix from a nested list.
 108 fromList :: (Vector v (w a), Vector w a, Vector v a) => [[a]] -> Mat v w a
 109 fromList vs = Mat (V.fromList $ map V.fromList vs)
 110
 111
 112 -- | Unsafe indexing.
 113 (!!!) :: (Vector w a) => Mat v w a -> (Int, Int) -> a
 114 (!!!) m (i, j) = (row m i) ! j
 115
 116 -- | Safe indexing.
 117 (!!?) :: Mat v w a -> (Int, Int) -> Maybe a
 118 (!!?) m@(Mat rows) (i, j)
 119   | i < 0 || j < 0 = Nothing
 120   | i > V.length rows = Nothing
 121   | otherwise = if j > V.length (row m j)
 122                 then Nothing
 123                 else Just $ (row m j) ! j
 124
 125
 126 -- | The number of rows in the matrix.
 127 nrows :: Mat v w a -> Int
 128 nrows (Mat rows) = V.length rows
 129
 130 -- | The number of columns in the first row of the
 131 --   matrix. Implementation stolen from Data.Vector.Fixed.length.
 132 ncols :: forall v w a. (Vector w a) => Mat v w a -> Int
 133 ncols _ = (arity (undefined :: Dim w))
 134
 135 -- | Return the @i@th row of @m@. Unsafe.
 136 row :: Mat v w a -> Int -> w a
 137 row (Mat rows) i = rows ! i
 138
 139
 140 -- | Return the @j@th column of @m@. Unsafe.
 141 column :: (Vector v a) => Mat v w a -> Int -> v a
 142 column (Mat rows) j =
 143   V.map (element j) rows
 144   where
 145     element = flip (!)
 146
 147
 148 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
 149 --   implementation.. it would be a little cleaner to use imap, but it
 150 --   doesn't seem to work.
 151 --
 152 --   TODO: Don't cheat with fromList.
 153 --
 154 --   Examples:
 155 --
 156 --   >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
 157 --   >>> transpose m
 158 --   ((1,3),(2,4))
 159 --
 160 transpose :: (Vector w (v a),
 161               Vector v a,
 162               Vector w a)
 163              => Mat v w a
 164              -> Mat w v a
 165 transpose m = Mat $ V.fromList column_list
 166   where
 167     column_list = [ column m i | i <- [0..(ncols m)-1] ]
 168
 169
 170 -- | Is @m@ symmetric?
 171 --
 172 --   Examples:
 173 --
 174 --   >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
 175 --   >>> symmetric m1
 176 --   True
 177 --
 178 --   >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
 179 --   >>> symmetric m2
 180 --   False
 181 --
 182 symmetric :: (Vector v (w a),
 183               Vector w a,
 184               v ~ w,
 185               Vector w Bool,
 186               Eq a)
 187              => Mat v w a
 188              -> Bool
 189 symmetric m =
 190   m == (transpose m)
 191
 192
 193 -- | Construct a new matrix from a function @lambda@. The function
 194 --   @lambda@ should take two parameters i,j corresponding to the
 195 --   entries in the matrix. The i,j entry of the resulting matrix will
 196 --   have the value returned by lambda i j.
 197 --
 198 --   TODO: Don't cheat with fromList.
 199 --
 200 --   Examples:
 201 --
 202 --   >>> let lambda i j = i + j
 203 --   >>> construct lambda :: Mat3 Int
 204 --   ((0,1,2),(1,2,3),(2,3,4))
 205 --
 206 construct :: forall v w a.
 207              (Vector v (w a),
 208               Vector w a)
 209              => (Int -> Int -> a)
 210              -> Mat v w a
 211 construct lambda = Mat rows
 212   where
 213     -- The arity trick is used in Data.Vector.Fixed.length.
 214     imax = (arity (undefined :: Dim v)) - 1
 215     jmax = (arity (undefined :: Dim w)) - 1
 216     row' i = V.fromList [ lambda i j | j <- [0..jmax] ]
 217     rows = V.fromList [ row' i | i <- [0..imax] ]
 218
 219 -- | Given a positive-definite matrix @m@, computes the
 220 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 221 --   values on the diagonal of @r@ positive.
 222 --
 223 --   Examples:
 224 --
 225 --   >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
 226 --   >>> cholesky m1
 227 --   ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))
 228 --   >>> (transpose (cholesky m1)) `mult` (cholesky m1)
 229 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 230 --
 231 cholesky :: forall a v w.
 232             (Algebraic.C a,
 233              Vector v (w a),
 234              Vector w a,
 235              Vector v a)
 236             => (Mat v w a)
 237             -> (Mat v w a)
 238 cholesky m = construct r
 239   where
 240     r :: Int -> Int -> a
 241     r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
 242           | i < j =
 243               (((m !!! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i)
 244           | otherwise = 0
 245
 246
 247 -- | Matrix multiplication. Our 'Num' instance doesn't define one, and
 248 --   we need additional restrictions on the result type anyway.
 249 --
 250 --   Examples:
 251 --
 252 --   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat Vec2D Vec3D Int
 253 --   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat Vec3D Vec2D Int
 254 --   >>> m1 `mult` m2
 255 --   ((22,28),(49,64))
 256 --
 257 mult :: (Ring.C a,
 258          Vector v a,
 259          Vector w a,
 260          Vector z a,
 261          Vector v (z a))
 262         => Mat v w a
 263         -> Mat w z a
 264         -> Mat v z a
 265 mult m1 m2 = construct lambda
 266   where
 267     lambda i j =
 268       sum [(m1 !!! (i,k)) * (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
 269
 270
 271
 272 instance (Ring.C a,
 273           Vector v (w a),
 274           Vector w a)
 275          => Additive.C (Mat v w a) where
 276
 277   (Mat rows1) + (Mat rows2) =
 278     Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
 279
 280   (Mat rows1) - (Mat rows2) =
 281     Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
 282
 283   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 284
 285
 286 instance (Ring.C a,
 287           Vector v (w a),
 288           Vector w a,
 289           v ~ w)
 290          => Ring.C (Mat v w a) where
 291     one = Mat (V.replicate $ V.replicate (fromInteger 1))
 292     m1 * m2 = m1 `mult` m2
 293
 294
 295 instance (Algebraic.C a,
 296           ToRational.C a,
 297           Vector v (w a),
 298           Vector w a,
 299           Vector v a,
 300           Vector v [a])
 301          => Normed (Mat v w a) where
 302    -- Treat the matrix as a big vector.
 303    norm_p p (Mat rows) =
 304      sqrt $ sum [(fromRational' $ toRational x)^2 | x <- xs]
 305      where
 306        xs = concat $ V.toList $ V.map V.toList rows
 307
 308    norm_infty m@(Mat rows)
 309      | nrows m == 0 || ncols m == 0 = 0
 310      | otherwise =
 311          fromRational' $ toRational $
 312          P.maximum $ V.toList $ V.map (P.maximum . V.toList) rows
 313
 314
 315
 316
 317
 318 -- Vector helpers. We want it to be easy to create low-dimension
 319 -- column vectors.
 320 type Vec a b = Mat a D1 b
 321
 322 vec2d :: (a,a) -> Mat D2 D1 a
 323 vec2d (x,y) = Mat (D2 (D1 x) (D1 y))
 324
 325 vec3d :: (a,a,a) -> Mat D3 D1 a
 326 vec3d (x,y,z) = Mat (D3 (D1 x) (D1 y) (D1 z))
 327
 328 vec4d :: (a,a,a,a) -> Mat D4 D1 a
 329 vec4d (w,x,y,z) = Mat (D4 (D1 w) (D1 x) (D1 y) (D1 z))
 330
 331 dot :: (RealRing.C a,
 332         Dim w ~ V.N1,
 333         Vector v a,
 334         Vector w a,
 335         Vector w (v a),
 336         Vector w (w a))
 337        => Mat v w a
 338        -> Mat v w a
 339        -> a
 340 v1 `dot` v2 = ((transpose v1) `mult` v2) !!! (0, 0)
 341
 342
 343 -- | The angle between @v1@ and @v2@ in Euclidean space.
 344 --
 345 --   Examples:
 346 --
 347 --   >>> let v1 = make2d (1.0, 0.0)
 348 --   >>> let v2 = make2d (0.0, 1.0)
 349 --   >>> angle v1 v2 == pi/2.0
 350 --   True
 351 --
 352 angle :: (Transcendental.C a,
 353           RealRing.C a,
 354           Dim w ~ V.N1,
 355           Vector w (w a),
 356           Vector v [a],
 357           Vector v a,
 358           Vector w a,
 359           Vector v (w a),
 360           Vector w (v a),
 361           ToRational.C a)
 362           => Mat v w a
 363           -> Mat v w a
 364           -> a
 365 angle v1 v2 =
 366   acos theta
 367   where
 368    theta = (recip norms) * (v1 `dot` v2)
 369    norms = (norm v1) * (norm v2)