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