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