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   generate,
  23   mk0,
  24   mk1,
  25   mk2,
  26   mk3,
  27   mk4,
  28   mk5 )
  29 import qualified Data.Vector.Fixed as V (
  30   and,
  31   foldl,
  32   fromList,
  33   head,
  34   ifoldl,
  35   ifoldr,
  36   imap,
  37   map,
  38   replicate,
  39   reverse,
  40   toList,
  41   zipWith )
  42 import Data.Vector.Fixed.Cont ( Arity, arity )
  43 import Linear.Vector ( Vec, delete )
  44 import Naturals
  45 import Normed ( Normed(..) )
  46
  47 -- We want the "max" that works on Ord, not the one that only works on
  48 -- Bool/Integer from the Lattice class!
  49 import NumericPrelude hiding ( (*), abs, max)
  50 import qualified NumericPrelude as NP ( (*) )
  51 import qualified Algebra.Absolute as Absolute ( C )
  52 import Algebra.Absolute ( abs )
  53 import qualified Algebra.Additive as Additive ( C )
  54 import qualified Algebra.Algebraic as Algebraic ( C )
  55 import Algebra.Algebraic ( root )
  56 import qualified Algebra.Field as Field ( C )
  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, max)
  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 Mat0 a = Mat Z Z a
  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 type Mat6 a = Mat N6 N6 a
  75 type Mat7 a = Mat N7 N7 a
  76
  77 -- * Type synonyms for 1-by-n row "vectors".
  78
  79 -- | Type synonym for row vectors expressed as 1-by-n matrices.
  80 type Row n a = Mat N1 n a
  81
  82 type Row1 a = Row N1 a
  83 type Row2 a = Row N2 a
  84 type Row3 a = Row N3 a
  85 type Row4 a = Row N4 a
  86 type Row5 a = Row N5 a
  87
  88 -- * Type synonyms for n-by-1 column "vectors".
  89
  90 -- | Type synonym for column vectors expressed as n-by-1 matrices.
  91 type Col n a = Mat n N1 a
  92
  93 type Col0 a = Col Z a
  94 type Col1 a = Col N1 a
  95 type Col2 a = Col N2 a
  96 type Col3 a = Col N3 a
  97 type Col4 a = Col N4 a
  98 type Col5 a = Col N5 a
  99 type Col6 a = Col N6 a
 100 type Col7 a = Col N7 a
 101 type Col8 a = Col N8 a
 102 type Col9 a = Col N9 a
 103 type Col10 a = Col N10 a
 104 type Col11 a = Col N11 a
 105 type Col12 a = Col N12 a
 106 type Col13 a = Col N13 a
 107 type Col14 a = Col N14 a
 108 type Col15 a = Col N15 a
 109 type Col16 a = Col N16 a
 110 type Col17 a = Col N17 a
 111 type Col18 a = Col N18 a
 112 type Col19 a = Col N19 a
 113 type Col20 a = Col N20 a
 114 type Col21 a = Col N21 a
 115 type Col22 a = Col N22 a
 116 type Col23 a = Col N23 a
 117 type Col24 a = Col N24 a
 118 type Col25 a = Col N25 a
 119 type Col26 a = Col N26 a
 120 type Col27 a = Col N27 a
 121 type Col28 a = Col N28 a
 122 type Col29 a = Col N29 a
 123 type Col30 a = Col N30 a
 124 type Col31 a = Col N31 a
 125 type Col32 a = Col N32 a
 126
 127
 128 instance (Eq a) => Eq (Mat m n a) where
 129   -- | Compare a row at a time.
 130   --
 131   --   Examples:
 132   --
 133   --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
 134   --   >>> let m2 = fromList [[1,2],[3,4]] :: Mat2 Int
 135   --   >>> let m3 = fromList [[5,6],[7,8]] :: Mat2 Int
 136   --   >>> m1 == m2
 137   --   True
 138   --   >>> m1 == m3
 139   --   False
 140   --
 141   (Mat rows_one) == (Mat rows_two) =
 142     V.and $ V.zipWith comp rows_one rows_two
 143     where
 144       -- Compare a row, one column at a time.
 145       comp row1 row2 = V.and (V.zipWith (==) row1 row2)
 146
 147
 148 instance (Show a) => Show (Mat m n a) where
 149   -- | Display matrices and vectors as ordinary tuples. This is poor
 150   --   practice, but these results are primarily displayed
 151   --   interactively and convenience trumps correctness (said the guy
 152   --   who insists his vector lengths be statically checked at
 153   --   compile-time).
 154   --
 155   --   Examples:
 156   --
 157   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 158   --   >>> show m
 159   --   ((1,2),(3,4))
 160   --
 161   show (Mat rows) =
 162     "(" ++ (intercalate "," (V.toList row_strings)) ++ ")"
 163     where
 164       row_strings = V.map show_vector rows
 165       show_vector v1 =
 166         "(" ++ (intercalate "," element_strings) ++ ")"
 167         where
 168           v1l = V.toList v1
 169           element_strings = P.map show v1l
 170
 171
 172 -- | Convert a matrix to a nested list.
 173 toList :: Mat m n a -> [[a]]
 174 toList (Mat rows) = map V.toList (V.toList rows)
 175
 176
 177 -- | Create a matrix from a nested list.
 178 fromList :: (Arity m, Arity n) => [[a]] -> Mat m n a
 179 fromList vs = Mat (V.fromList $ map V.fromList vs)
 180
 181
 182 -- | Unsafe indexing. Much faster than the safe indexing.
 183 (!!!) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> a
 184 (!!!) (Mat rows) (i, j) = (rows ! i) ! j
 185
 186
 187 -- | Safe indexing.
 188 --
 189 --   Examples:
 190 --
 191 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 192 --   >>> m !!? (-1,-1)
 193 --   Nothing
 194 --   >>>  m !!? (-1,0)
 195 --   Nothing
 196 --   >>>  m !!? (-1,1)
 197 --   Nothing
 198 --   >>>  m !!? (0,-1)
 199 --   Nothing
 200 --   >>>  m !!? (0,0)
 201 --   Just 1
 202 --   >>>  m !!? (0,1)
 203 --   Just 2
 204 --   >>>  m !!? (1,-1)
 205 --   Nothing
 206 --   >>>  m !!? (1,0)
 207 --   Just 3
 208 --   >>>  m !!? (1,1)
 209 --   Just 4
 210 --   >>>  m !!? (2,-1)
 211 --   Nothing
 212 --   >>>  m !!? (2,0)
 213 --   Nothing
 214 --   >>>  m !!? (2,1)
 215 --   Nothing
 216 --   >>>  m !!? (2,2)
 217 --   Nothing
 218 --
 219 (!!?) :: (Arity m, Arity n) => Mat m n a -> (Int, Int) -> Maybe a
 220 (!!?) matrix idx =
 221   ifoldl2 f Nothing matrix
 222   where
 223     f k l found cur = if (k,l) == idx then (Just cur) else found
 224
 225
 226 -- | The number of rows in the matrix.
 227 nrows :: forall m n a. (Arity m) => Mat m n a -> Int
 228 nrows _ = arity (undefined :: m)
 229
 230
 231 -- | The number of columns in the first row of the
 232 --   matrix. Implementation stolen from Data.Vector.Fixed.length.
 233 ncols :: forall m n a. (Arity n) => Mat m n a -> Int
 234 ncols _ = arity (undefined :: n)
 235
 236
 237 -- | Return the @i@th row of @m@ as a matrix. Unsafe.
 238 row :: (Arity m, Arity n) => Mat m n a -> Int -> Row n a
 239 row m i =
 240   construct lambda
 241   where
 242     lambda _ j = m !!! (i, j)
 243
 244
 245 -- | Return the @j@th column of @m@ as a matrix. Unsafe.
 246 column :: (Arity m, Arity n) => Mat m n a -> Int -> Col m a
 247 column m j =
 248   construct lambda
 249   where
 250     lambda i _ = m !!! (i, j)
 251
 252
 253 -- | Transpose @m@; switch it's columns and its rows. This is a dirty
 254 --   implementation, but I don't see a better way.
 255 --
 256 --   TODO: Don't cheat with fromList.
 257 --
 258 --   Examples:
 259 --
 260 --   >>> let m = fromList [[1,2], [3,4]] :: Mat2 Int
 261 --   >>> transpose m
 262 --   ((1,3),(2,4))
 263 --
 264 transpose :: (Arity m, Arity n) => Mat m n a -> Mat n m a
 265 transpose matrix =
 266   construct lambda
 267   where
 268     lambda i j = matrix !!! (j,i)
 269
 270
 271 -- | Is @m@ symmetric?
 272 --
 273 --   Examples:
 274 --
 275 --   >>> let m1 = fromList [[1,2], [2,1]] :: Mat2 Int
 276 --   >>> symmetric m1
 277 --   True
 278 --
 279 --   >>> let m2 = fromList [[1,2], [3,1]] :: Mat2 Int
 280 --   >>> symmetric m2
 281 --   False
 282 --
 283 symmetric :: (Eq a, Arity m) => Mat m m a -> Bool
 284 symmetric m =
 285   m == (transpose m)
 286
 287
 288 -- | Construct a new matrix from a function @lambda@. The function
 289 --   @lambda@ should take two parameters i,j corresponding to the
 290 --   entries in the matrix. The i,j entry of the resulting matrix will
 291 --   have the value returned by lambda i j.
 292 --
 293 --   Examples:
 294 --
 295 --   >>> let lambda i j = i + j
 296 --   >>> construct lambda :: Mat3 Int
 297 --   ((0,1,2),(1,2,3),(2,3,4))
 298 --
 299 construct :: forall m n a. (Arity m, Arity n)
 300           => (Int -> Int -> a) -> Mat m n a
 301 construct lambda = Mat $ generate make_row
 302   where
 303     make_row :: Int -> Vec n a
 304     make_row i = generate (lambda i)
 305
 306
 307 -- | Create an identity matrix with the right dimensions.
 308 --
 309 --   Examples:
 310 --
 311 --   >>> identity_matrix :: Mat3 Int
 312 --   ((1,0,0),(0,1,0),(0,0,1))
 313 --   >>> identity_matrix :: Mat3 Double
 314 --   ((1.0,0.0,0.0),(0.0,1.0,0.0),(0.0,0.0,1.0))
 315 --
 316 identity_matrix :: (Arity m, Ring.C a) => Mat m m a
 317 identity_matrix =
 318   construct (\i j -> if i == j then (fromInteger 1) else (fromInteger 0))
 319
 320
 321 -- | Given a positive-definite matrix @m@, computes the
 322 --   upper-triangular matrix @r@ with (transpose r)*r == m and all
 323 --   values on the diagonal of @r@ positive.
 324 --
 325 --   Examples:
 326 --
 327 --   >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
 328 --   >>> let r = cholesky m1
 329 --   >>> frobenius_norm ((transpose r)*r - m1) < 1e-10
 330 --   True
 331 --   >>> is_upper_triangular r
 332 --   True
 333 --
 334 --   >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
 335 --   >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
 336 --   >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
 337 --   >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
 338 --   >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
 339 --   >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
 340 --   >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
 341 --   >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
 342 --
 343 --   >>> let e1 = [2.449489742783178,0,0,0,0,0,0] :: [Double]
 344 --   >>> let e2 = [-1.224744871391589,3,0,0,0,0,0] :: [Double]
 345 --   >>> let e3 = [0,-5/2,5/2,0,0,0,0] :: [Double]
 346 --   >>> let e4 = [0,0,0,2.449489742783178,0,0,0] :: [Double]
 347 --   >>> let e5 = [0,0,0,0,2.449489742783178,0,0] :: [Double]
 348 --   >>> let e6 = [0,0,0,0,0,2.449489742783178,0] :: [Double]
 349 --   >>> let e7 = [0,0,0,0,0,0,3.872983346207417] :: [Double]
 350 --   >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Mat N7 N7 Double
 351 --
 352 --   >>> let r = cholesky big_K
 353 --   >>> frobenius_norm (r - (transpose expected)) < 1e-12
 354 --   True
 355 --
 356 cholesky :: forall m a. (Algebraic.C a, Arity m)
 357          => (Mat m m a) -> (Mat m m a)
 358 cholesky m = ifoldl2 f zero m
 359   where
 360     f :: Int -> Int -> (Mat m m a) -> a -> (Mat m m a)
 361     f i j cur_R _ = set_idx cur_R (i,j) (r cur_R i j)
 362
 363     r :: (Mat m m a) -> Int -> Int -> a
 364     r cur_R i j
 365       | i == j = sqrt(m !!! (i,j) - sum [(cur_R !!! (k,i))^2 | k <- [0..i-1]])
 366       | i < j = (((m !!! (i,j))
 367                   - sum [(cur_R !!! (k,i)) NP.* (cur_R !!! (k,j))
 368                              | k <- [0..i-1]]))/(cur_R !!! (i,i))
 369       | otherwise = 0
 370
 371
 372
 373 -- | Returns True if the given matrix is upper-triangular, and False
 374 --   otherwise. The parameter @epsilon@ lets the caller choose a
 375 --   tolerance.
 376 --
 377 --   Examples:
 378 --
 379 --   >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
 380 --   >>> is_upper_triangular m
 381 --   False
 382 --   >>> is_upper_triangular' 1e-10 m
 383 --   True
 384 --
 385 is_upper_triangular' :: forall m n a.
 386                         (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
 387                     => a -- ^ The tolerance @epsilon@.
 388                     -> Mat m n a
 389                     -> Bool
 390 is_upper_triangular' epsilon matrix =
 391   ifoldl2 f True matrix
 392   where
 393     f :: Int -> Int -> Bool -> a -> Bool
 394     f _ _ False _ = False
 395     f i j True x
 396       | i <= j = True
 397       -- use "less than or equal to" so zero is a valid epsilon
 398       | otherwise = abs x <= epsilon
 399
 400
 401 -- | Returns True if the given matrix is upper-triangular, and False
 402 --   otherwise. We don't delegate to the general
 403 --   'is_upper_triangular'' here because it imposes additional
 404 --   typeclass constraints throughout the library.
 405 --
 406 --   Examples:
 407 --
 408 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 409 --   >>> is_upper_triangular m
 410 --   False
 411 --
 412 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 413 --   >>> is_upper_triangular m
 414 --   True
 415 --
 416 is_upper_triangular :: forall m n a.
 417                        (Eq a, Ring.C a, Arity m, Arity n)
 418                     => Mat m n a -> Bool
 419 is_upper_triangular matrix =
 420   ifoldl2 f True matrix
 421   where
 422     f :: Int -> Int -> Bool -> a -> Bool
 423     f _ _ False _ = False
 424     f i j True x
 425       | i <= j = True
 426       | otherwise = x == 0
 427
 428
 429
 430 -- | Returns True if the given matrix is lower-triangular, and False
 431 --   otherwise.
 432 --
 433 --   Examples:
 434 --
 435 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 436 --   >>> is_lower_triangular m
 437 --   True
 438 --
 439 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 440 --   >>> is_lower_triangular m
 441 --   False
 442 --
 443 is_lower_triangular :: (Eq a,
 444                         Ring.C a,
 445                         Arity m,
 446                         Arity n)
 447                     => Mat m n a
 448                     -> Bool
 449 is_lower_triangular = is_upper_triangular . transpose
 450
 451
 452 -- | Returns True if the given matrix is lower-triangular, and False
 453 --   otherwise. The parameter @epsilon@ lets the caller choose a
 454 --   tolerance.
 455 --
 456 --   Examples:
 457 --
 458 --   >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
 459 --   >>> is_lower_triangular m
 460 --   False
 461 --   >>> is_lower_triangular' 1e-12 m
 462 --   True
 463 --
 464 is_lower_triangular' :: (Ord a,
 465                          Ring.C a,
 466                          Absolute.C a,
 467                          Arity m,
 468                          Arity n)
 469                     => a -- ^ The tolerance @epsilon@.
 470                     -> Mat m n a
 471                     -> Bool
 472 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
 473
 474
 475 -- | Returns True if the given matrix is triangular, and False
 476 --   otherwise.
 477 --
 478 --   Examples:
 479 --
 480 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 481 --   >>> is_triangular m
 482 --   True
 483 --
 484 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 485 --   >>> is_triangular m
 486 --   True
 487 --
 488 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 489 --   >>> is_triangular m
 490 --   False
 491 --
 492 is_triangular :: (Ord a,
 493                   Ring.C a,
 494                   Absolute.C a,
 495                   Arity m,
 496                   Arity n)
 497               => Mat m n a
 498               -> Bool
 499 is_triangular m = is_upper_triangular m || is_lower_triangular m
 500
 501
 502 -- | Delete the @i@th row and @j@th column from the matrix. The name
 503 --   \"preminor\" is made up, but is meant to signify that this is
 504 --   usually used in the computationof a minor. A minor is simply the
 505 --   determinant of a preminor in that case.
 506 --
 507 --   Examples:
 508 --
 509 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 510 --   >>> preminor m 0 0 :: Mat2 Int
 511 --   ((5,6),(8,9))
 512 --   >>> preminor m 1 1 :: Mat2 Int
 513 --   ((1,3),(7,9))
 514 --
 515 preminor :: (Arity m, Arity n)
 516       => Mat (S m) (S n) a
 517       -> Int
 518       -> Int
 519       -> Mat m n a
 520 preminor (Mat rows) i j = m
 521   where
 522     rows' = delete rows i
 523     m = Mat $ V.map ((flip delete) j) rows'
 524
 525
 526 -- | Compute the i,jth minor of a @matrix@.
 527 --
 528 --   Examples:
 529 --
 530 --   >>> let m1 = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Double
 531 --   >>> minor m1 1 1
 532 --   -12.0
 533 --
 534 minor :: (Arity m, Determined (Mat m m) a)
 535       => Mat (S m) (S m) a
 536       -> Int
 537       -> Int
 538       -> a
 539 minor matrix i j = determinant (preminor matrix i j)
 540
 541 class (Eq a, Ring.C a) => Determined p a where
 542   determinant :: (p a) -> a
 543
 544 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
 545   determinant = unscalar
 546
 547 instance (Ord a,
 548           Ring.C a,
 549           Absolute.C a,
 550           Arity n,
 551           Determined (Mat (S n) (S n)) a)
 552          => Determined (Mat (S (S n)) (S (S n))) a where
 553   -- | The recursive definition with a special-case for triangular matrices.
 554   --
 555   --   Examples:
 556   --
 557   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 558   --   >>> determinant m
 559   --   -1
 560   --
 561   determinant m
 562     | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
 563     | otherwise = determinant_recursive
 564         where
 565           m' i j = m !!! (i,j)
 566
 567           determinant_recursive =
 568             sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (minor m 0 j)
 569               | j <- [0..(ncols m)-1] ]
 570
 571
 572
 573 -- | Matrix multiplication.
 574 --
 575 --   Examples:
 576 --
 577 --   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat N2 N3 Int
 578 --   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
 579 --   >>> m1 * m2
 580 --   ((22,28),(49,64))
 581 --
 582 infixl 7 *
 583 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
 584         => Mat (S m) (S n) a
 585         -> Mat (S n) (S p) a
 586         -> Mat (S m) (S p) a
 587 (*) m1 m2 = construct lambda
 588   where
 589     lambda i j = (transpose $ row m1 i) `dot` (column m2 j)
 590
 591
 592
 593 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
 594
 595   (Mat rows_one) + (Mat rows_two) =
 596     Mat $ V.zipWith (V.zipWith (+)) rows_one rows_two
 597
 598   (Mat rows_one) - (Mat rows_two) =
 599     Mat $ V.zipWith (V.zipWith (-)) rows_one rows_two
 600
 601   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 602
 603
 604 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat (S m) (S n) a) where
 605   -- The first * is ring multiplication, the second is matrix
 606   -- multiplication.
 607   one = identity_matrix
 608   m1 * m2 = m1 * m2
 609
 610
 611 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
 612   -- We can multiply a matrix by a scalar of the same type as its
 613   -- elements.
 614   x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
 615
 616
 617 instance (Absolute.C a,
 618           Algebraic.C a,
 619           ToRational.C a,
 620           Arity m)
 621          => Normed (Col m a) where
 622   -- | Generic p-norms for vectors in R^m that are represented as m-by-1
 623   --   matrices.
 624   --
 625   --   Examples:
 626   --
 627   --   >>> let v1 = vec2d (3,4)
 628   --   >>> norm_p 1 v1
 629   --   7.0
 630   --   >>> norm_p 2 v1
 631   --   5.0
 632   --
 633   --   >>> let v1 = vec2d (-1,1) :: Col2 Double
 634   --   >>> norm_p 1 v1 :: Double
 635   --   2.0
 636   --
 637   --   >>> let v1 = vec0d :: Col0 Double
 638   --   >>> norm v1
 639   --   0.0
 640   --
 641   norm_p p (Mat rows) =
 642     (root p') $ sum [fromRational' (toRational $ abs x)^p' | x <- xs]
 643     where
 644       p' = toInteger p
 645       xs = concat $ V.toList $ V.map V.toList rows
 646
 647   -- | The infinity norm.
 648   --
 649   --   Examples:
 650   --
 651   --   >>> let v1 = vec3d (1,5,2)
 652   --   >>> norm_infty v1
 653   --   5
 654   --
 655   norm_infty (Mat rows) =
 656     fromRational' $ toRational
 657                       $ (V.foldl P.max 0) $ V.map (V.foldl P.max 0) rows
 658
 659
 660 -- | Compute the Frobenius norm of a matrix. This essentially treats
 661 --   the matrix as one long vector containing all of its entries (in
 662 --   any order, it doesn't matter).
 663 --
 664 --   Examples:
 665 --
 666 --   >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
 667 --   >>> frobenius_norm m == sqrt 285
 668 --   True
 669 --
 670 --   >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
 671 --   >>> frobenius_norm m == 3
 672 --   True
 673 --
 674 frobenius_norm :: (Arity m, Arity n, Algebraic.C a, Ring.C a)
 675                => Mat m n a
 676                -> a
 677 frobenius_norm matrix =
 678   sqrt $ element_sum2 $ squares
 679   where
 680     squares = map2 (^2) matrix
 681
 682
 683 -- Vector helpers. We want it to be easy to create low-dimension
 684 -- column vectors, which are nx1 matrices.
 685
 686 -- | Convenient constructor for 2D vectors.
 687 --
 688 --   Examples:
 689 --
 690 --   >>> import Roots.Simple
 691 --   >>> let fst m = m !!! (0,0)
 692 --   >>> let snd m = m !!! (1,0)
 693 --   >>> let h = 0.5 :: Double
 694 --   >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
 695 --   >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
 696 --   >>> let g u = vec2d ((g1 u), (g2 u))
 697 --   >>> let u0 = vec2d (1.0, 1.0)
 698 --   >>> let eps = 1/(10^9)
 699 --   >>> fixed_point g eps u0
 700 --   ((1.0728549599342185),(1.0820591495686167))
 701 --
 702 vec0d :: Col0 a
 703 vec0d = Mat mk0
 704
 705 vec1d :: (a) -> Col1 a
 706 vec1d (x) = Mat (mk1 (mk1 x))
 707
 708 vec2d :: (a,a) -> Col2 a
 709 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
 710
 711 vec3d :: (a,a,a) -> Col3 a
 712 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
 713
 714 vec4d :: (a,a,a,a) -> Col4 a
 715 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 716
 717 vec5d :: (a,a,a,a,a) -> Col5 a
 718 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 719
 720
 721 -- Since we commandeered multiplication, we need to create 1x1
 722 -- matrices in order to multiply things.
 723 scalar :: a -> Mat1 a
 724 scalar x = Mat (mk1 (mk1 x))
 725
 726 -- Get the scalar value out of a 1x1 matrix.
 727 unscalar :: Mat1 a -> a
 728 unscalar (Mat rows) = V.head $ V.head rows
 729
 730
 731 dot :: (Ring.C a, Arity m)
 732        => Col (S m) a
 733        -> Col (S m) a
 734        -> a
 735 v1 `dot` v2 = element_sum2 $ zipwith2 (NP.*) v1 v2
 736
 737
 738 -- | The angle between @v1@ and @v2@ in Euclidean space.
 739 --
 740 --   Examples:
 741 --
 742 --   >>> let v1 = vec2d (1.0, 0.0)
 743 --   >>> let v2 = vec2d (0.0, 1.0)
 744 --   >>> angle v1 v2 == pi/2.0
 745 --   True
 746 --
 747 angle :: (Transcendental.C a,
 748           RealRing.C a,
 749           m ~ S t,
 750           Arity t,
 751           ToRational.C a)
 752           => Col m a
 753           -> Col m a
 754           -> a
 755 angle v1 v2 =
 756   acos theta
 757   where
 758    theta = (recip norms) NP.* (v1 `dot` v2)
 759    norms = (norm v1) NP.* (norm v2)
 760
 761
 762 -- | Retrieve the diagonal elements of the given matrix as a \"column
 763 --   vector,\" i.e. a m-by-1 matrix. We require the matrix to be
 764 --   square to avoid ambiguity in the return type which would ideally
 765 --   have dimension min(m,n) supposing an m-by-n matrix.
 766 --
 767 --   Examples:
 768 --
 769 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 770 --   >>> diagonal m
 771 --   ((1),(5),(9))
 772 --
 773 diagonal :: (Arity m) => Mat m m a -> Col m a
 774 diagonal matrix =
 775   construct lambda
 776   where
 777     lambda i _ = matrix !!! (i,i)
 778
 779
 780 -- | Given a square @matrix@, return a new matrix of the same size
 781 --   containing only the on-diagonal entries of @matrix@. The
 782 --   off-diagonal entries are set to zero.
 783 --
 784 --   Examples:
 785 --
 786 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 787 --   >>> diagonal_part m
 788 --   ((1,0,0),(0,5,0),(0,0,9))
 789 --
 790 diagonal_part :: (Arity m, Ring.C a)
 791          => Mat m m a
 792          -> Mat m m a
 793 diagonal_part matrix =
 794   construct lambda
 795   where
 796     lambda i j = if i == j then matrix !!! (i,j) else 0
 797
 798
 799 -- | Given a square @matrix@, return a new matrix of the same size
 800 --   containing only the on-diagonal and below-diagonal entries of
 801 --   @matrix@. The above-diagonal entries are set to zero.
 802 --
 803 --   Examples:
 804 --
 805 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 806 --   >>> lt_part m
 807 --   ((1,0,0),(4,5,0),(7,8,9))
 808 --
 809 lt_part :: (Arity m, Ring.C a)
 810         => Mat m m a
 811         -> Mat m m a
 812 lt_part matrix =
 813   construct lambda
 814   where
 815     lambda i j = if i >= j then matrix !!! (i,j) else 0
 816
 817
 818 -- | Given a square @matrix@, return a new matrix of the same size
 819 --   containing only the below-diagonal entries of @matrix@. The on-
 820 --   and above-diagonal entries are set to zero.
 821 --
 822 --   Examples:
 823 --
 824 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 825 --   >>> lt_part_strict m
 826 --   ((0,0,0),(4,0,0),(7,8,0))
 827 --
 828 lt_part_strict :: (Arity m, Ring.C a)
 829         => Mat m m a
 830         -> Mat m m a
 831 lt_part_strict matrix =
 832   construct lambda
 833   where
 834     lambda i j = if i > j then matrix !!! (i,j) else 0
 835
 836
 837 -- | Given a square @matrix@, return a new matrix of the same size
 838 --   containing only the on-diagonal and above-diagonal entries of
 839 --   @matrix@. The below-diagonal entries are set to zero.
 840 --
 841 --   Examples:
 842 --
 843 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 844 --   >>> ut_part m
 845 --   ((1,2,3),(0,5,6),(0,0,9))
 846 --
 847 ut_part :: (Arity m, Ring.C a)
 848         => Mat m m a
 849         -> Mat m m a
 850 ut_part = transpose . lt_part . transpose
 851
 852
 853 -- | Given a square @matrix@, return a new matrix of the same size
 854 --   containing only the above-diagonal entries of @matrix@. The on-
 855 --   and below-diagonal entries are set to zero.
 856 --
 857 --   Examples:
 858 --
 859 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 860 --   >>> ut_part_strict m
 861 --   ((0,2,3),(0,0,6),(0,0,0))
 862 --
 863 ut_part_strict :: (Arity m, Ring.C a)
 864         => Mat m m a
 865         -> Mat m m a
 866 ut_part_strict = transpose . lt_part_strict . transpose
 867
 868
 869 -- | Compute the trace of a square matrix, the sum of the elements
 870 --   which lie on its diagonal. We require the matrix to be
 871 --   square to avoid ambiguity in the return type which would ideally
 872 --   have dimension min(m,n) supposing an m-by-n matrix.
 873 --
 874 --   Examples:
 875 --
 876 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 877 --   >>> trace m
 878 --   15
 879 --
 880 trace :: (Arity m, Ring.C a) => Mat m m a -> a
 881 trace = element_sum2 . diagonal
 882
 883
 884
 885 -- | Zip together two matrices.
 886 --
 887 --   TODO: don't cheat with construct (map V.zips instead).
 888 --
 889 --   Examples:
 890 --
 891 --   >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
 892 --   >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
 893 --   >>> zip2 m1 m2
 894 --   (((1,1)),((1,2)),((1,3)))
 895 --
 896 --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
 897 --   >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
 898 --   >>> zip2 m1 m2
 899 --   (((1,1),(2,1)),((3,1),(4,1)))
 900 --
 901 zip2 :: (Arity m, Arity n) => Mat m n a -> Mat m n b -> Mat m n (a,b)
 902 zip2 m1 m2 =
 903   construct lambda
 904   where
 905     lambda i j = (m1 !!! (i,j), m2 !!! (i,j))
 906
 907
 908 -- | Zip together three matrices.
 909 --
 910 --   TODO: don't cheat with construct (map V.zips instead).
 911 --
 912 --   Examples:
 913 --
 914 --   >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
 915 --   >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
 916 --   >>> let m3 = fromList [[4],[5],[6]] :: Col3 Int
 917 --   >>> zip2three m1 m2 m3
 918 --   (((1,1,4)),((1,2,5)),((1,3,6)))
 919 --
 920 --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
 921 --   >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
 922 --   >>> let m3 = fromList [[8,2],[6,3]] :: Mat2 Int
 923 --   >>> zip2three m1 m2 m3
 924 --   (((1,1,8),(2,1,2)),((3,1,6),(4,1,3)))
 925 --
 926 zip2three :: (Arity m, Arity n)
 927           => Mat m n a
 928           -> Mat m n a
 929           -> Mat m n a
 930           -> Mat m n (a,a,a)
 931 zip2three m1 m2 m3 =
 932   construct lambda
 933   where
 934     lambda i j = (m1 !!! (i,j), m2 !!! (i,j), m3 !!! (i,j))
 935
 936
 937 -- | Zip together two matrices using the supplied function.
 938 --
 939 --   Examples:
 940 --
 941 --   >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
 942 --   >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
 943 --   >>> zipwith2 (^) c1 c2
 944 --   ((1),(32),(729))
 945 --
 946 zipwith2 :: (Arity m, Arity n)
 947            => (a -> b -> c)
 948            -> Mat m n a
 949            -> Mat m n b
 950            -> Mat m n c
 951 zipwith2 f c1 c2 =
 952   construct lambda
 953   where
 954     lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
 955
 956
 957 -- | Map a function over a matrix of any dimensions.
 958 --
 959 --   Examples:
 960 --
 961 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 962 --   >>> map2 (^2) m
 963 --   ((1,4),(9,16))
 964 --
 965 map2 :: (a -> b) -> Mat m n a -> Mat m n b
 966 map2 f (Mat rows) =
 967   Mat $ V.map g rows
 968   where
 969     g = V.map f
 970
 971
 972 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
 973 --   (of the row/column) to the accumulation function. The fold occurs
 974 --   from top-left to bottom-right.
 975 --
 976 --   Examples:
 977 --
 978 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 979 --   >>> ifoldl2 (\i j cur _ -> cur + i + j) 0 m
 980 --   18
 981 --
 982 ifoldl2 :: forall a b m n.
 983            (Int -> Int -> b -> a -> b)
 984         -> b
 985         -> Mat m n a
 986         -> b
 987 ifoldl2 f initial (Mat rows) =
 988   V.ifoldl row_function initial rows
 989   where
 990     -- | The order that we need this in (so that @g idx@ makes sense)
 991     --   is a little funny. So that we don't need to pass weird
 992     --   functions into ifoldl2, we swap the second and third
 993     --   arguments of @f@ calling the result @g@.
 994     g :: Int -> b -> Int -> a -> b
 995     g w x y = f w y x
 996
 997     row_function :: b -> Int -> Vec n a -> b
 998     row_function rowinit idx r = V.ifoldl (g idx) rowinit r
 999
1000
1001 -- | Left fold over the entries of a matrix (top-left to bottom-right).
1002 --
1003 foldl2 :: forall a b m n.
1004           (b -> a -> b)
1005         -> b
1006         -> Mat m n a
1007         -> b
1008 foldl2 f initial matrix =
1009   -- Use the index fold but ignore the index arguments.
1010   let g _ _ = f in ifoldl2 g initial matrix
1011
1012
1013 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
1014 --   (of the row/column) to the accumulation function. The fold occurs
1015 --   from bottom-right to top-left.
1016 --
1017 --   The order of the arguments in the supplied function are different
1018 --   from those in V.ifoldr; we keep them similar to ifoldl2.
1019 --
1020 --   Examples:
1021 --
1022 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1023 --   >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
1024 --   18
1025 --
1026 ifoldr2 :: forall a b m n.
1027            (Int -> Int -> b -> a -> b)
1028         -> b
1029         -> Mat m n a
1030         -> b
1031 ifoldr2 f initial (Mat rows) =
1032   V.ifoldr row_function initial rows
1033   where
1034     -- | Swap the order of arguments in @f@ so that it agrees with the
1035     --   @f@ passed to ifoldl2.
1036     g :: Int -> Int -> a -> b -> b
1037     g w x y z = f w x z y
1038
1039     row_function :: Int -> Vec n a -> b -> b
1040     row_function idx r rowinit = V.ifoldr (g idx) rowinit r
1041
1042
1043 -- | Map a function over a matrix of any dimensions, passing the
1044 --   coordinates @i@ and @j@ to the function @f@.
1045 --
1046 --   Examples:
1047 --
1048 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
1049 --   >>> imap2 (\i j _ -> i+j) m
1050 --   ((0,1),(1,2))
1051 --
1052 imap2 :: (Int -> Int -> a -> b) -> Mat m n a -> Mat m n b
1053 imap2 f (Mat rows) =
1054   Mat $ V.imap g rows
1055   where
1056     g i = V.imap (f i)
1057
1058
1059 -- | Reverse the order of elements in a matrix.
1060 --
1061 --   Examples:
1062 --
1063 --   >>> let m1 = fromList [[1,2,3]] :: Row3 Int
1064 --   >>> reverse2 m1
1065 --   ((3,2,1))
1066 --
1067 --   >>> let m1 = vec3d (1,2,3 :: Int)
1068 --   >>> reverse2 m1
1069 --   ((3),(2),(1))
1070 --
1071 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1072 --   >>> reverse2 m
1073 --   ((9,8,7),(6,5,4),(3,2,1))
1074 --
1075 reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
1076 reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
1077
1078
1079 -- | Unsafely set the (i,j) element of the given matrix.
1080 --
1081 --   Examples:
1082 --
1083 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1084 --   >>> set_idx m (1,1) 17
1085 --   ((1,2,3),(4,17,6),(7,8,9))
1086 --
1087 set_idx :: forall m n a.
1088            (Arity m, Arity n)
1089         => Mat m n a
1090         -> (Int, Int)
1091         -> a
1092         -> Mat m n a
1093 set_idx matrix (i,j) newval =
1094   imap2 updater matrix
1095   where
1096     updater :: Int -> Int -> a -> a
1097     updater k l existing =
1098       if k == i && l == j
1099       then newval
1100       else existing
1101
1102
1103 -- | Compute the i,jth cofactor of the given @matrix@. This simply
1104 --   premultiplues the i,jth minor by (-1)^(i+j).
1105 cofactor :: (Arity m, Determined (Mat m m) a)
1106          => Mat (S m) (S m) a
1107          -> Int
1108          -> Int
1109          -> a
1110 cofactor matrix i j =
1111   (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j)
1112
1113
1114 -- | Compute the inverse of a matrix using cofactor expansion
1115 --   (generalized Cramer's rule).
1116 --
1117 --   Examples:
1118 --
1119 --   >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double
1120 --   >>> let e1 = [54/1624, -22/1624] :: [Double]
1121 --   >>> let e2 = [-17/1624, 37/1624] :: [Double]
1122 --   >>> let expected = fromList [e1, e2] :: Mat2 Double
1123 --   >>> let actual = inverse m1
1124 --   >>> frobenius_norm (actual - expected) < 1e-12
1125 --   True
1126 --
1127 inverse :: (Arity m,
1128             Determined (Mat (S m) (S m)) a,
1129             Determined (Mat m m) a,
1130             Field.C a)
1131         => Mat (S m) (S m) a
1132         -> Mat (S m) (S m) a
1133 inverse matrix =
1134   (1 / (determinant matrix)) *> (transpose $ construct lambda)
1135   where
1136     lambda i j = cofactor matrix i j
1137
1138
1139
1140 -- | Retrieve the rows of a matrix as a column matrix. If the given
1141 --   matrix is m-by-n, the result would be an m-by-1 column whose
1142 --   entries are 1-by-n row matrices.
1143 --
1144 --   Examples:
1145 --
1146 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
1147 --   >>> (rows2 m) !!! (0,0)
1148 --   ((1,2))
1149 --   >>> (rows2 m) !!! (1,0)
1150 --   ((3,4))
1151 --
1152 rows2 :: (Arity m, Arity n)
1153       => Mat m n a
1154       -> Col m (Row n a)
1155 rows2 (Mat rows) =
1156   Mat $ V.map (mk1. Mat . mk1) rows
1157
1158
1159
1160 -- | Sum the elements of a matrix.
1161 --
1162 --   Examples:
1163 --
1164 --   >>> let m = fromList [[1,-1],[3,4]] :: Mat2 Int
1165 --   >>> element_sum2 m
1166 --   7
1167 --
1168 element_sum2 :: (Arity m, Arity n, Additive.C a) => Mat m n a -> a
1169 element_sum2 = foldl2 (+) zero