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