src/Linear/Matrix.hs

   1 {-# LANGUAGE ExistentialQuantification #-}
   2 {-# LANGUAGE FlexibleContexts #-}
   3 {-# LANGUAGE FlexibleInstances #-}
   4 {-# LANGUAGE MultiParamTypeClasses #-}
   5 {-# LANGUAGE NoMonomorphismRestriction #-}
   6 {-# LANGUAGE ScopedTypeVariables #-}
   7 {-# LANGUAGE TypeFamilies #-}
   8 {-# LANGUAGE RebindableSyntax #-}
   9
  10 -- | Boxed matrices; that is, boxed m-vectors of boxed n-vectors. We
  11 --   assume that the underlying representation is
  12 --   Data.Vector.Fixed.Boxed.Vec for simplicity. It was tried in
  13 --   generality and failed.
  14 --
  15 module Linear.Matrix
  16 where
  17
  18 import Data.List (intercalate)
  19
  20 import Data.Vector.Fixed (
  21   (!),
  22   N1,
  23   N2,
  24   N3,
  25   N4,
  26   N5,
  27   S,
  28   Z,
  29   generate,
  30   mk1,
  31   mk2,
  32   mk3,
  33   mk4,
  34   mk5
  35   )
  36 import qualified Data.Vector.Fixed as V (
  37   and,
  38   fromList,
  39   head,
  40   length,
  41   map,
  42   maximum,
  43   replicate,
  44   toList,
  45   zipWith )
  46 import Data.Vector.Fixed.Cont ( Arity, arity )
  47 import Linear.Vector ( Vec, delete, element_sum )
  48 import Normed ( Normed(..) )
  49
  50 import NumericPrelude hiding ( (*), abs )
  51 import qualified NumericPrelude as NP ( (*) )
  52 import qualified Algebra.Absolute as Absolute ( C )
  53 import Algebra.Absolute ( abs )
  54 import qualified Algebra.Additive as Additive ( C )
  55 import qualified Algebra.Algebraic as Algebraic ( C )
  56 import Algebra.Algebraic ( root )
  57 import qualified Algebra.Ring as Ring ( C )
  58 import qualified Algebra.Module as Module ( C )
  59 import qualified Algebra.RealRing as RealRing ( C )
  60 import qualified Algebra.ToRational as ToRational ( C )
  61 import qualified Algebra.Transcendental as Transcendental ( C )
  62 import qualified Prelude as P ( map )
  63
  64 -- | Our main matrix type.
  65 data Mat m n a = (Arity m, Arity n) => Mat (Vec m (Vec n a))
  66
  67 -- Type synonyms for n-by-n matrices.
  68 type Mat1 a = Mat N1 N1 a
  69 type Mat2 a = Mat N2 N2 a
  70 type Mat3 a = Mat N3 N3 a
  71 type Mat4 a = Mat N4 N4 a
  72 type Mat5 a = Mat N5 N5 a
  73
  74 -- | Type synonym for row vectors expressed as 1-by-n matrices.
  75 type Row n a = Mat N1 n a
  76
  77 -- Type synonyms for 1-by-n row "vectors".
  78 type Row1 a = Row N1 a
  79 type Row2 a = Row N2 a
  80 type Row3 a = Row N3 a
  81 type Row4 a = Row N4 a
  82 type Row5 a = Row N5 a
  83
  84 -- | Type synonym for column vectors expressed as n-by-1 matrices.
  85 type Col n a = Mat n N1 a
  86
  87 -- Type synonyms for n-by-1 column "vectors".
  88 type Col1 a = Col N1 a
  89 type Col2 a = Col N2 a
  90 type Col3 a = Col N3 a
  91 type Col4 a = Col N4 a
  92 type Col5 a = Col N5 a
  93
  94 instance (Eq a) => Eq (Mat m n a) where
  95   -- | Compare a row at a time.
  96   --
  97   --   Examples:
  98   --
  99   --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
 100   --   >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
 101   --   >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
 102   --   >>> m1 == m2
 103   --   True
 104   --   >>> m1 == m3
 105   --   False
 106   --
 107   (Mat rows1) == (Mat rows2) =
 108     V.and $ V.zipWith comp rows1 rows2
 109     where
 110       -- Compare a row, one column at a time.
 111       comp row1 row2 = V.and (V.zipWith (==) row1 row2)
 112
 113
 114 instance (Show a) => Show (Mat m n a) where
 115   -- | Display matrices and vectors as ordinary tuples. This is poor
 116   --   practice, but these results are primarily displayed
 117   --   interactively and convenience trumps correctness (said the guy
 118   --   who insists his vector lengths be statically checked at
 119   --   compile-time).
 120   --
 121   --   Examples:
 122   --
 123   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 124   --   >>> show m
 125   --   ((1,2),(3,4))
 126   --
 127   show (Mat rows) =
 128     "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
 129     where
 130       row_strings = V.map show_vector rows
 131       show_vector v1 =
 132         "(" ++ (intercalate "," element_strings) ++ ")"
 133         where
 134           v1l = V.toList v1
 135           element_strings = P.map show v1l
 136
 137
 138 -- | Convert a matrix to a nested list.
 139 toList :: Mat m n a -> [[a]]
 140 toList (Mat rows) = map V.toList (V.toList rows)
 141
 142 -- | Create a matrix from a nested list.
 143 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
 144 fromList vs = Mat (V.fromList $ map V.fromList vs)
 145
 146
 147 -- | Unsafe indexing.
 148 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
 149 (!!!) m (i, j) = (row m i) ! j
 150
 151 -- | Safe indexing.
 152 (!!?) :: Mat m n a -> (Int, Int) -> Maybe a
 153 (!!?) m@(Mat rows) (i, j)
 154   | i < 0 || j < 0 = Nothing
 155   | i > V.length rows = Nothing
 156   | otherwise = if j > V.length (row m j)
 157                 then Nothing
 158                 else Just $ (row m j) ! j
 159
 160
 161 -- | The number of rows in the matrix.
 162 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
 163 nrows _ = arity (undefined :: m)
 164
 165 -- | The number of columns in the first row of the
 166 --   matrix. Implementation stolen from Data.Vector.Fixed.length.
 167 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
 168 ncols _ = arity (undefined :: n)
 169
 170
 171 -- | Return the @i@th row of @m@. Unsafe.
 172 row :: Mat m n a -> Int -> (Vec n a)
 173 row (Mat rows) i = rows ! i
 174
 175
 176 -- | Return the @i@th row of @m@ as a matrix. Unsafe.
 177 row' :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
 178 row' m i =
 179   construct lambda
 180   where
 181     lambda _ j = m !!! (i, j)
 182
 183
 184 -- | Return the @j@th column of @m@. Unsafe.
 185 column :: Mat m n a -> Int -> (Vec m a)
 186 column (Mat rows) j =
 187   V.map (element j) rows
 188   where
 189     element = flip (!)
 190
 191
 192 -- | Return the @j@th column of @m@ as a matrix. Unsafe.
 193 column' :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
 194 column' m j =
 195   construct lambda
 196   where
 197     lambda i _ = m !!! (i, j)
 198
 199
 200 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
 201 --   implementation.. it would be a little cleaner to use imap, but it
 202 --   doesn't seem to work.
 203 --
 204 --   TODO: Don't cheat with fromList.
 205 --
 206 --   Examples:
 207 --
 208 --   >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
 209 --   >>> transpose m
 210 --   ((1,3),(2,4))
 211 --
 212 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
 213 transpose m = Mat $ V.fromList column_list
 214   where
 215     column_list = [ column m i | i <- [0..(ncols m)-1] ]
 216
 217
 218 -- | Is @m@ symmetric?
 219 --
 220 --   Examples:
 221 --
 222 --   >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
 223 --   >>> symmetric m1
 224 --   True
 225 --
 226 --   >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
 227 --   >>> symmetric m2
 228 --   False
 229 --
 230 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
 231 symmetric m =
 232   m == (transpose m)
 233
 234
 235 -- | Construct a new matrix from a function @lambda@. The function
 236 --   @lambda@ should take two parameters i,j corresponding to the
 237 --   entries in the matrix. The i,j entry of the resulting matrix will
 238 --   have the value returned by lambda i j.
 239 --
 240 --   Examples:
 241 --
 242 --   >>> let lambda i j = i + j
 243 --   >>> construct lambda :: Mat3 Int
 244 --   ((0,1,2),(1,2,3),(2,3,4))
 245 --
 246 construct :: forall m n a. (Arity m, Arity n)
 247           => (Int -> Int -> a) -> Mat m n a
 248 construct lambda = Mat $ generate make_row
 249   where
 250     make_row :: Int -> Vec n a
 251     make_row i = generate (lambda i)
 252
 253
 254 -- | Create an identity matrix with the right dimensions.
 255 --
 256 --   Examples:
 257 --
 258 --   >>> identity_matrix :: Mat3 Int
 259 --   ((1,0,0),(0,1,0),(0,0,1))
 260 --   >>> identity_matrix :: Mat3 Double
 261 --   ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
 262 --
 263 identity_matrix :: (Arity m, Ring.C a) => Mat m m a
 264 identity_matrix =
 265   construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
 266
 267 -- | Given a positive-definite matrix @m@, computes the
 268 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 269 --   values on the diagonal of @r@ positive.
 270 --
 271 --   Examples:
 272 --
 273 --   >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
 274 --   >>> cholesky m1
 275 --   ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))
 276 --   >>> (transpose (cholesky m1)) * (cholesky m1)
 277 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 278 --
 279 cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
 280          => (Mat m n a) -> (Mat m n a)
 281 cholesky m = construct r
 282   where
 283     r :: Int -> Int -> a
 284     r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
 285           | i < j =
 286               (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
 287           | otherwise = 0
 288
 289
 290 -- | Returns True if the given matrix is upper-triangular, and False
 291 --   otherwise. The parameter @epsilon@ lets the caller choose a
 292 --   tolerance.
 293 --
 294 --   Examples:
 295 --
 296 --   >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
 297 --   >>> is_upper_triangular m
 298 --   False
 299 --   >>> is_upper_triangular' 1e-10 m
 300 --   True
 301 --
 302 --   TODO:
 303 --
 304 --     1. Don't cheat with lists.
 305 --
 306 is_upper_triangular' :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
 307                     => a -- ^ The tolerance @epsilon@.
 308                     -> Mat m n a
 309                     -> Bool
 310 is_upper_triangular' epsilon m =
 311   and $ concat results
 312   where
 313     results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
 314
 315     test :: Int -> Int -> Bool
 316     test i j
 317       | i <= j = True
 318       -- use "less than or equal to" so zero is a valid epsilon
 319       | otherwise = abs (m !!! (i,j)) <= epsilon
 320
 321
 322 -- | Returns True if the given matrix is upper-triangular, and False
 323 --   otherwise. A specialized version of 'is_upper_triangular\'' with
 324 --   @epsilon = 0@.
 325 --
 326 --   Examples:
 327 --
 328 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 329 --   >>> is_upper_triangular m
 330 --   False
 331 --
 332 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 333 --   >>> is_upper_triangular m
 334 --   True
 335 --
 336 --   TODO:
 337 --
 338 --     1. The Ord constraint is too strong here, Eq would suffice.
 339 --
 340 is_upper_triangular :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
 341                     => Mat m n a -> Bool
 342 is_upper_triangular = is_upper_triangular' 0
 343
 344
 345 -- | Returns True if the given matrix is lower-triangular, and False
 346 --   otherwise. This is a specialized version of 'is_lower_triangular\''
 347 --   with @epsilon = 0@.
 348 --
 349 --   Examples:
 350 --
 351 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 352 --   >>> is_lower_triangular m
 353 --   True
 354 --
 355 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 356 --   >>> is_lower_triangular m
 357 --   False
 358 --
 359 is_lower_triangular :: (Ord a,
 360                         Ring.C a,
 361                         Absolute.C a,
 362                         Arity m,
 363                         Arity n)
 364                     => Mat m n a
 365                     -> Bool
 366 is_lower_triangular = is_upper_triangular . transpose
 367
 368
 369 -- | Returns True if the given matrix is lower-triangular, and False
 370 --   otherwise. The parameter @epsilon@ lets the caller choose a
 371 --   tolerance.
 372 --
 373 --   Examples:
 374 --
 375 --   >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
 376 --   >>> is_lower_triangular m
 377 --   False
 378 --   >>> is_lower_triangular' 1e-12 m
 379 --   True
 380 --
 381 is_lower_triangular' :: (Ord a,
 382                          Ring.C a,
 383                          Absolute.C a,
 384                          Arity m,
 385                          Arity n)
 386                     => a -- ^ The tolerance @epsilon@.
 387                     -> Mat m n a
 388                     -> Bool
 389 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
 390
 391
 392 -- | Returns True if the given matrix is triangular, and False
 393 --   otherwise.
 394 --
 395 --   Examples:
 396 --
 397 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 398 --   >>> is_triangular m
 399 --   True
 400 --
 401 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 402 --   >>> is_triangular m
 403 --   True
 404 --
 405 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 406 --   >>> is_triangular m
 407 --   False
 408 --
 409 is_triangular :: (Ord a,
 410                   Ring.C a,
 411                   Absolute.C a,
 412                   Arity m,
 413                   Arity n)
 414               => Mat m n a
 415               -> Bool
 416 is_triangular m = is_upper_triangular m || is_lower_triangular m
 417
 418
 419 -- | Return the (i,j)th minor of m.
 420 --
 421 --   Examples:
 422 --
 423 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 424 --   >>> minor m 0 0 :: Mat2 Int
 425 --   ((5,6),(8,9))
 426 --   >>> minor m 1 1 :: Mat2 Int
 427 --   ((1,3),(7,9))
 428 --
 429 minor :: (m ~ S r,
 430           n ~ S t,
 431           Arity r,
 432           Arity t)
 433       => Mat m n a
 434       -> Int
 435       -> Int
 436       -> Mat r t a
 437 minor (Mat rows) i j = m
 438   where
 439     rows' = delete rows i
 440     m = Mat $ V.map ((flip delete) j) rows'
 441
 442
 443 class (Eq a, Ring.C a) => Determined p a where
 444   determinant :: (p a) -> a
 445
 446 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
 447   determinant (Mat rows) = (V.head . V.head) rows
 448
 449 instance (Ord a,
 450           Ring.C a,
 451           Absolute.C a,
 452           Arity n,
 453           Determined (Mat (S n) (S n)) a)
 454          => Determined (Mat (S (S n)) (S (S n))) a where
 455   -- | The recursive definition with a special-case for triangular matrices.
 456   --
 457   --   Examples:
 458   --
 459   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 460   --   >>> determinant m
 461   --   -1
 462   --
 463   determinant m
 464     | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
 465     | otherwise = determinant_recursive
 466         where
 467           m' i j = m !!! (i,j)
 468
 469           det_minor i j = determinant (minor m i j)
 470
 471           determinant_recursive =
 472             sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j)
 473               | j <- [0..(ncols m)-1] ]
 474
 475
 476
 477 -- | Matrix multiplication.
 478 --
 479 --   Examples:
 480 --
 481 --   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat N2 N3 Int
 482 --   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
 483 --   >>> m1 * m2
 484 --   ((22,28),(49,64))
 485 --
 486 infixl 7 *
 487 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
 488         => Mat m n a
 489         -> Mat n p a
 490         -> Mat m p a
 491 (*) m1 m2 = construct lambda
 492   where
 493     lambda i j =
 494       sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
 495
 496
 497
 498 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
 499
 500   (Mat rows1) + (Mat rows2) =
 501     Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
 502
 503   (Mat rows1) - (Mat rows2) =
 504     Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
 505
 506   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 507
 508
 509 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where
 510   -- The first * is ring multiplication, the second is matrix
 511   -- multiplication.
 512   m1 * m2 = m1 * m2
 513
 514
 515 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
 516   -- We can multiply a matrix by a scalar of the same type as its
 517   -- elements.
 518   x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
 519
 520
 521 instance (Algebraic.C a,
 522           ToRational.C a,
 523           Arity m)
 524          => Normed (Mat (S m) N1 a) where
 525   -- | Generic p-norms for vectors in R^n that are represented as nx1
 526   --   matrices.
 527   --
 528   --   Examples:
 529   --
 530   --   >>> let v1 = vec2d (3,4)
 531   --   >>> norm_p 1 v1
 532   --   7.0
 533   --   >>> norm_p 2 v1
 534   --   5.0
 535   --
 536   norm_p p (Mat rows) =
 537     (root p') $ sum [fromRational' (toRational x)^p' | x <- xs]
 538     where
 539       p' = toInteger p
 540       xs = concat $ V.toList $ V.map V.toList rows
 541
 542   -- | The infinity norm.
 543   --
 544   --   Examples:
 545   --
 546   --   >>> let v1 = vec3d (1,5,2)
 547   --   >>> norm_infty v1
 548   --   5
 549   --
 550   norm_infty (Mat rows) =
 551     fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
 552
 553
 554 -- | Compute the Frobenius norm of a matrix. This essentially treats
 555 --   the matrix as one long vector containing all of its entries (in
 556 --   any order, it doesn't matter).
 557 --
 558 --   Examples:
 559 --
 560 --   >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
 561 --   >>> frobenius_norm m == sqrt 285
 562 --   True
 563 --
 564 --   >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
 565 --   >>> frobenius_norm m == 3
 566 --   True
 567 --
 568 frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
 569 frobenius_norm (Mat rows) =
 570   sqrt $ element_sum $ V.map row_sum rows
 571   where
 572     -- | Square and add up the entries of a row.
 573     row_sum = element_sum . V.map (^2)
 574
 575
 576 -- Vector helpers. We want it to be easy to create low-dimension
 577 -- column vectors, which are nx1 matrices.
 578
 579 -- | Convenient constructor for 2D vectors.
 580 --
 581 --   Examples:
 582 --
 583 --   >>> import Roots.Simple
 584 --   >>> let fst m = m !!! (0,0)
 585 --   >>> let snd m = m !!! (1,0)
 586 --   >>> let h = 0.5 :: Double
 587 --   >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
 588 --   >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
 589 --   >>> let g u = vec2d ((g1 u), (g2 u))
 590 --   >>> let u0 = vec2d (1.0, 1.0)
 591 --   >>> let eps = 1/(10^9)
 592 --   >>> fixed_point g eps u0
 593 --   ((1.0728549599342185),(1.0820591495686167))
 594 --
 595 vec1d :: (a) -> Col1 a
 596 vec1d (x) = Mat (mk1 (mk1 x))
 597
 598 vec2d :: (a,a) -> Col2 a
 599 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
 600
 601 vec3d :: (a,a,a) -> Col3 a
 602 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
 603
 604 vec4d :: (a,a,a,a) -> Col4 a
 605 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 606
 607 vec5d :: (a,a,a,a,a) -> Col5 a
 608 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 609
 610 -- Since we commandeered multiplication, we need to create 1x1
 611 -- matrices in order to multiply things.
 612 scalar :: a -> Mat1 a
 613 scalar x = Mat (mk1 (mk1 x))
 614
 615 dot :: (RealRing.C a, n ~ N1, m ~ S t, Arity t)
 616        => Mat m n a
 617        -> Mat m n a
 618        -> a
 619 v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)
 620
 621
 622 -- | The angle between @v1@ and @v2@ in Euclidean space.
 623 --
 624 --   Examples:
 625 --
 626 --   >>> let v1 = vec2d (1.0, 0.0)
 627 --   >>> let v2 = vec2d (0.0, 1.0)
 628 --   >>> angle v1 v2 == pi/2.0
 629 --   True
 630 --
 631 angle :: (Transcendental.C a,
 632           RealRing.C a,
 633           n ~ N1,
 634           m ~ S t,
 635           Arity t,
 636           ToRational.C a)
 637           => Mat m n a
 638           -> Mat m n a
 639           -> a
 640 angle v1 v2 =
 641   acos theta
 642   where
 643    theta = (recip norms) NP.* (v1 `dot` v2)
 644    norms = (norm v1) NP.* (norm v2)
 645
 646
 647 -- | Retrieve the diagonal elements of the given matrix as a \"column
 648 --   vector,\" i.e. a m-by-1 matrix. We require the matrix to be
 649 --   square to avoid ambiguity in the return type which would ideally
 650 --   have dimension min(m,n) supposing an m-by-n matrix.
 651 --
 652 --   Examples:
 653 --
 654 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 655 --   >>> diagonal m
 656 --   ((1),(5),(9))
 657 --
 658 diagonal :: (Arity m) => Mat m m a -> Col m a
 659 diagonal matrix =
 660   construct lambda
 661   where
 662     lambda i _ = matrix !!! (i,i)
 663
 664
 665 -- | Given a square @matrix@, return a new matrix of the same size
 666 --   containing only the on-diagonal entries of @matrix@. The
 667 --   off-diagonal entries are set to zero.
 668 --
 669 --   Examples:
 670 --
 671 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 672 --   >>> diagonal_part m
 673 --   ((1,0,0),(0,5,0),(0,0,9))
 674 --
 675 diagonal_part :: (Arity m, Ring.C a)
 676          => Mat m m a
 677          -> Mat m m a
 678 diagonal_part matrix =
 679   construct lambda
 680   where
 681     lambda i j = if i == j then matrix !!! (i,j) else 0
 682
 683
 684 -- | Given a square @matrix@, return a new matrix of the same size
 685 --   containing only the on-diagonal and below-diagonal entries of
 686 --   @matrix@. The above-diagonal entries are set to zero.
 687 --
 688 --   Examples:
 689 --
 690 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 691 --   >>> lt_part m
 692 --   ((1,0,0),(4,5,0),(7,8,9))
 693 --
 694 lt_part :: (Arity m, Ring.C a)
 695         => Mat m m a
 696         -> Mat m m a
 697 lt_part matrix =
 698   construct lambda
 699   where
 700     lambda i j = if i >= j then matrix !!! (i,j) else 0
 701
 702
 703 -- | Given a square @matrix@, return a new matrix of the same size
 704 --   containing only the below-diagonal entries of @matrix@. The on-
 705 --   and above-diagonal entries are set to zero.
 706 --
 707 --   Examples:
 708 --
 709 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 710 --   >>> lt_part_strict m
 711 --   ((0,0,0),(4,0,0),(7,8,0))
 712 --
 713 lt_part_strict :: (Arity m, Ring.C a)
 714         => Mat m m a
 715         -> Mat m m a
 716 lt_part_strict matrix =
 717   construct lambda
 718   where
 719     lambda i j = if i > j then matrix !!! (i,j) else 0
 720
 721
 722 -- | Given a square @matrix@, return a new matrix of the same size
 723 --   containing only the on-diagonal and above-diagonal entries of
 724 --   @matrix@. The below-diagonal entries are set to zero.
 725 --
 726 --   Examples:
 727 --
 728 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 729 --   >>> ut_part m
 730 --   ((1,2,3),(0,5,6),(0,0,9))
 731 --
 732 ut_part :: (Arity m, Ring.C a)
 733         => Mat m m a
 734         -> Mat m m a
 735 ut_part = transpose . lt_part . transpose
 736
 737
 738 -- | Given a square @matrix@, return a new matrix of the same size
 739 --   containing only the above-diagonal entries of @matrix@. The on-
 740 --   and below-diagonal entries are set to zero.
 741 --
 742 --   Examples:
 743 --
 744 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 745 --   >>> ut_part_strict m
 746 --   ((0,2,3),(0,0,6),(0,0,0))
 747 --
 748 ut_part_strict :: (Arity m, Ring.C a)
 749         => Mat m m a
 750         -> Mat m m a
 751 ut_part_strict = transpose . lt_part_strict . transpose
 752
 753
 754 -- | Compute the trace of a square matrix, the sum of the elements
 755 --   which lie on its diagonal. We require the matrix to be
 756 --   square to avoid ambiguity in the return type which would ideally
 757 --   have dimension min(m,n) supposing an m-by-n matrix.
 758 --
 759 --   Examples:
 760 --
 761 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 762 --   >>> trace m
 763 --   15
 764 --
 765 trace :: (Arity m, Ring.C a) => Mat m m a -> a
 766 trace matrix =
 767   let (Mat rows) = diagonal matrix
 768   in
 769     element_sum $ V.map V.head rows