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