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, element_sum )
  43 import Naturals ( N1, N2, N3, N4, N5, N6, N7, N8, N9, N10, S, Z )
  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
  71 -- * Type synonyms for 1-by-n row "vectors".
  72
  73 -- | Type synonym for row vectors expressed as 1-by-n matrices.
  74 type Row n a = Mat N1 n a
  75
  76 type Row1 a = Row N1 a
  77 type Row2 a = Row N2 a
  78 type Row3 a = Row N3 a
  79 type Row4 a = Row N4 a
  80 type Row5 a = Row N5 a
  81
  82 -- * Type synonyms for n-by-1 column "vectors".
  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 Col1 a = Col N1 a
  88 type Col2 a = Col N2 a
  89 type Col3 a = Col N3 a
  90 type Col4 a = Col N4 a
  91 type Col5 a = Col N5 a
  92 type Col6 a = Col N6 a
  93 type Col7 a = Col N7 a
  94 type Col8 a = Col N8 a
  95 type Col9 a = Col N9 a
  96 type Col10 a = Col N10 a -- We need a big column for Gaussian quadrature.
  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 --   >>> let r = cholesky m1
 300 --   >>> frobenius_norm ((transpose r)*r - m1) < 1e-10
 301 --   True
 302 --   >>> is_upper_triangular r
 303 --   True
 304 --
 305 --   >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
 306 --   >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
 307 --   >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
 308 --   >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
 309 --   >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
 310 --   >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
 311 --   >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
 312 --   >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
 313 --
 314 --   >>> let e1 = [2.449489742783178,0,0,0,0,0,0] :: [Double]
 315 --   >>> let e2 = [-1.224744871391589,3,0,0,0,0,0] :: [Double]
 316 --   >>> let e3 = [0,-5/2,5/2,0,0,0,0] :: [Double]
 317 --   >>> let e4 = [0,0,0,2.449489742783178,0,0,0] :: [Double]
 318 --   >>> let e5 = [0,0,0,0,2.449489742783178,0,0] :: [Double]
 319 --   >>> let e6 = [0,0,0,0,0,2.449489742783178,0] :: [Double]
 320 --   >>> let e7 = [0,0,0,0,0,0,3.872983346207417] :: [Double]
 321 --   >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Mat N7 N7 Double
 322 --
 323 --   >>> let r = cholesky big_K
 324 --   >>> frobenius_norm (r - (transpose expected)) < 1e-12
 325 --   True
 326 --
 327 cholesky :: forall m n a. (Algebraic.C a, Arity m, Arity n)
 328          => (Mat m n a) -> (Mat m n a)
 329 cholesky m = construct r
 330   where
 331     r :: Int -> Int -> a
 332     r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
 333           | i < j =
 334               (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)
 335           | otherwise = 0
 336
 337
 338 -- | Returns True if the given matrix is upper-triangular, and False
 339 --   otherwise. The parameter @epsilon@ lets the caller choose a
 340 --   tolerance.
 341 --
 342 --   Examples:
 343 --
 344 --   >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
 345 --   >>> is_upper_triangular m
 346 --   False
 347 --   >>> is_upper_triangular' 1e-10 m
 348 --   True
 349 --
 350 is_upper_triangular' :: forall m n a.
 351                         (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
 352                     => a -- ^ The tolerance @epsilon@.
 353                     -> Mat m n a
 354                     -> Bool
 355 is_upper_triangular' epsilon matrix =
 356   ifoldl2 f True matrix
 357   where
 358     f :: Int -> Int -> Bool -> a -> Bool
 359     f _ _ False _ = False
 360     f i j True x
 361       | i <= j = True
 362       -- use "less than or equal to" so zero is a valid epsilon
 363       | otherwise = abs x <= epsilon
 364
 365
 366 -- | Returns True if the given matrix is upper-triangular, and False
 367 --   otherwise. We don't delegate to the general
 368 --   'is_upper_triangular'' here because it imposes additional
 369 --   typeclass constraints throughout the library.
 370 --
 371 --   Examples:
 372 --
 373 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 374 --   >>> is_upper_triangular m
 375 --   False
 376 --
 377 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 378 --   >>> is_upper_triangular m
 379 --   True
 380 --
 381 is_upper_triangular :: forall m n a.
 382                        (Eq a, Ring.C a, Arity m, Arity n)
 383                     => Mat m n a -> Bool
 384 is_upper_triangular matrix =
 385   ifoldl2 f True matrix
 386   where
 387     f :: Int -> Int -> Bool -> a -> Bool
 388     f _ _ False _ = False
 389     f i j True x
 390       | i <= j = True
 391       | otherwise = x == 0
 392
 393
 394
 395 -- | Returns True if the given matrix is lower-triangular, and False
 396 --   otherwise.
 397 --
 398 --   Examples:
 399 --
 400 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 401 --   >>> is_lower_triangular m
 402 --   True
 403 --
 404 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 405 --   >>> is_lower_triangular m
 406 --   False
 407 --
 408 is_lower_triangular :: (Eq a,
 409                         Ring.C a,
 410                         Arity m,
 411                         Arity n)
 412                     => Mat m n a
 413                     -> Bool
 414 is_lower_triangular = is_upper_triangular . transpose
 415
 416
 417 -- | Returns True if the given matrix is lower-triangular, and False
 418 --   otherwise. The parameter @epsilon@ lets the caller choose a
 419 --   tolerance.
 420 --
 421 --   Examples:
 422 --
 423 --   >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
 424 --   >>> is_lower_triangular m
 425 --   False
 426 --   >>> is_lower_triangular' 1e-12 m
 427 --   True
 428 --
 429 is_lower_triangular' :: (Ord a,
 430                          Ring.C a,
 431                          Absolute.C a,
 432                          Arity m,
 433                          Arity n)
 434                     => a -- ^ The tolerance @epsilon@.
 435                     -> Mat m n a
 436                     -> Bool
 437 is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
 438
 439
 440 -- | Returns True if the given matrix is triangular, and False
 441 --   otherwise.
 442 --
 443 --   Examples:
 444 --
 445 --   >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
 446 --   >>> is_triangular m
 447 --   True
 448 --
 449 --   >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
 450 --   >>> is_triangular m
 451 --   True
 452 --
 453 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 454 --   >>> is_triangular m
 455 --   False
 456 --
 457 is_triangular :: (Ord a,
 458                   Ring.C a,
 459                   Absolute.C a,
 460                   Arity m,
 461                   Arity n)
 462               => Mat m n a
 463               -> Bool
 464 is_triangular m = is_upper_triangular m || is_lower_triangular m
 465
 466
 467 -- | Delete the @i@th row and @j@th column from the matrix. The name
 468 --   \"preminor\" is made up, but is meant to signify that this is
 469 --   usually used in the computationof a minor. A minor is simply the
 470 --   determinant of a preminor in that case.
 471 --
 472 --   Examples:
 473 --
 474 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 475 --   >>> preminor m 0 0 :: Mat2 Int
 476 --   ((5,6),(8,9))
 477 --   >>> preminor m 1 1 :: Mat2 Int
 478 --   ((1,3),(7,9))
 479 --
 480 preminor :: (Arity m, Arity n)
 481       => Mat (S m) (S n) a
 482       -> Int
 483       -> Int
 484       -> Mat m n a
 485 preminor (Mat rows) i j = m
 486   where
 487     rows' = delete rows i
 488     m = Mat $ V.map ((flip delete) j) rows'
 489
 490
 491 -- | Compute the i,jth minor of a @matrix@.
 492 --
 493 --   Examples:
 494 --
 495 --   >>> let m1 = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Double
 496 --   >>> minor m1 1 1
 497 --   -12.0
 498 --
 499 minor :: (Arity m, Determined (Mat m m) a)
 500       => Mat (S m) (S m) a
 501       -> Int
 502       -> Int
 503       -> a
 504 minor matrix i j = determinant (preminor matrix i j)
 505
 506 class (Eq a, Ring.C a) => Determined p a where
 507   determinant :: (p a) -> a
 508
 509 instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
 510   determinant = unscalar
 511
 512 instance (Ord a,
 513           Ring.C a,
 514           Absolute.C a,
 515           Arity n,
 516           Determined (Mat (S n) (S n)) a)
 517          => Determined (Mat (S (S n)) (S (S n))) a where
 518   -- | The recursive definition with a special-case for triangular matrices.
 519   --
 520   --   Examples:
 521   --
 522   --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 523   --   >>> determinant m
 524   --   -1
 525   --
 526   determinant m
 527     | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
 528     | otherwise = determinant_recursive
 529         where
 530           m' i j = m !!! (i,j)
 531
 532           determinant_recursive =
 533             sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (minor m 0 j)
 534               | j <- [0..(ncols m)-1] ]
 535
 536
 537
 538 -- | Matrix multiplication.
 539 --
 540 --   Examples:
 541 --
 542 --   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat N2 N3 Int
 543 --   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat N3 N2 Int
 544 --   >>> m1 * m2
 545 --   ((22,28),(49,64))
 546 --
 547 infixl 7 *
 548 (*) :: (Ring.C a, Arity m, Arity n, Arity p)
 549         => Mat m n a
 550         -> Mat n p a
 551         -> Mat m p a
 552 (*) m1 m2 = construct lambda
 553   where
 554     lambda i j =
 555       sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
 556
 557
 558
 559 instance (Ring.C a, Arity m, Arity n) => Additive.C (Mat m n a) where
 560
 561   (Mat rows1) + (Mat rows2) =
 562     Mat $ V.zipWith (V.zipWith (+)) rows1 rows2
 563
 564   (Mat rows1) - (Mat rows2) =
 565     Mat $ V.zipWith (V.zipWith (-)) rows1 rows2
 566
 567   zero = Mat (V.replicate $ V.replicate (fromInteger 0))
 568
 569
 570 instance (Ring.C a, Arity m, Arity n, m ~ n) => Ring.C (Mat m n a) where
 571   -- The first * is ring multiplication, the second is matrix
 572   -- multiplication.
 573   m1 * m2 = m1 * m2
 574
 575
 576 instance (Ring.C a, Arity m, Arity n) => Module.C a (Mat m n a) where
 577   -- We can multiply a matrix by a scalar of the same type as its
 578   -- elements.
 579   x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows
 580
 581
 582 instance (Absolute.C a,
 583           Algebraic.C a,
 584           ToRational.C a,
 585           Arity m)
 586          => Normed (Col (S m) a) where
 587   -- | Generic p-norms for vectors in R^n that are represented as n-by-1
 588   --   matrices.
 589   --
 590   --   Examples:
 591   --
 592   --   >>> let v1 = vec2d (3,4)
 593   --   >>> norm_p 1 v1
 594   --   7.0
 595   --   >>> norm_p 2 v1
 596   --   5.0
 597   --
 598   --   >>> let v1 = vec2d (-1,1) :: Col2 Double
 599   --   >>> norm_p 1 v1 :: Double
 600   --   2.0
 601   --
 602   norm_p p (Mat rows) =
 603     (root p') $ sum [fromRational' (toRational $ abs x)^p' | x <- xs]
 604     where
 605       p' = toInteger p
 606       xs = concat $ V.toList $ V.map V.toList rows
 607
 608   -- | The infinity norm.
 609   --
 610   --   Examples:
 611   --
 612   --   >>> let v1 = vec3d (1,5,2)
 613   --   >>> norm_infty v1
 614   --   5
 615   --
 616   norm_infty (Mat rows) =
 617     fromRational' $ toRational $ V.maximum $ V.map V.maximum rows
 618
 619
 620 -- | Compute the Frobenius norm of a matrix. This essentially treats
 621 --   the matrix as one long vector containing all of its entries (in
 622 --   any order, it doesn't matter).
 623 --
 624 --   Examples:
 625 --
 626 --   >>> let m = fromList [[1, 2, 3],[4,5,6],[7,8,9]] :: Mat3 Double
 627 --   >>> frobenius_norm m == sqrt 285
 628 --   True
 629 --
 630 --   >>> let m = fromList [[1, -1, 1],[-1,1,-1],[1,-1,1]] :: Mat3 Double
 631 --   >>> frobenius_norm m == 3
 632 --   True
 633 --
 634 frobenius_norm :: (Algebraic.C a, Ring.C a) => Mat m n a -> a
 635 frobenius_norm (Mat rows) =
 636   sqrt $ element_sum $ V.map row_sum rows
 637   where
 638     -- | Square and add up the entries of a row.
 639     row_sum = element_sum . V.map (^2)
 640
 641
 642 -- Vector helpers. We want it to be easy to create low-dimension
 643 -- column vectors, which are nx1 matrices.
 644
 645 -- | Convenient constructor for 2D vectors.
 646 --
 647 --   Examples:
 648 --
 649 --   >>> import Roots.Simple
 650 --   >>> let fst m = m !!! (0,0)
 651 --   >>> let snd m = m !!! (1,0)
 652 --   >>> let h = 0.5 :: Double
 653 --   >>> let g1 m = 1.0 + h NP.* exp(-((fst m)^2))/(1.0 + (snd m)^2)
 654 --   >>> let g2 m = 0.5 + h NP.* atan((fst m)^2 + (snd m)^2)
 655 --   >>> let g u = vec2d ((g1 u), (g2 u))
 656 --   >>> let u0 = vec2d (1.0, 1.0)
 657 --   >>> let eps = 1/(10^9)
 658 --   >>> fixed_point g eps u0
 659 --   ((1.0728549599342185),(1.0820591495686167))
 660 --
 661 vec1d :: (a) -> Col1 a
 662 vec1d (x) = Mat (mk1 (mk1 x))
 663
 664 vec2d :: (a,a) -> Col2 a
 665 vec2d (x,y) = Mat (mk2 (mk1 x) (mk1 y))
 666
 667 vec3d :: (a,a,a) -> Col3 a
 668 vec3d (x,y,z) = Mat (mk3 (mk1 x) (mk1 y) (mk1 z))
 669
 670 vec4d :: (a,a,a,a) -> Col4 a
 671 vec4d (w,x,y,z) = Mat (mk4 (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 672
 673 vec5d :: (a,a,a,a,a) -> Col5 a
 674 vec5d (v,w,x,y,z) = Mat (mk5 (mk1 v) (mk1 w) (mk1 x) (mk1 y) (mk1 z))
 675
 676
 677 -- Since we commandeered multiplication, we need to create 1x1
 678 -- matrices in order to multiply things.
 679 scalar :: a -> Mat1 a
 680 scalar x = Mat (mk1 (mk1 x))
 681
 682 -- Get the scalar value out of a 1x1 matrix.
 683 unscalar :: Mat1 a -> a
 684 unscalar (Mat rows) = V.head $ V.head rows
 685
 686
 687 dot :: (Ring.C a, Arity m)
 688        => Col (S m) a
 689        -> Col (S m) a
 690        -> a
 691 v1 `dot` v2 = unscalar $ ((transpose v1) * v2)
 692
 693
 694 -- | The angle between @v1@ and @v2@ in Euclidean space.
 695 --
 696 --   Examples:
 697 --
 698 --   >>> let v1 = vec2d (1.0, 0.0)
 699 --   >>> let v2 = vec2d (0.0, 1.0)
 700 --   >>> angle v1 v2 == pi/2.0
 701 --   True
 702 --
 703 angle :: (Transcendental.C a,
 704           RealRing.C a,
 705           m ~ S t,
 706           Arity t,
 707           ToRational.C a)
 708           => Col m a
 709           -> Col m a
 710           -> a
 711 angle v1 v2 =
 712   acos theta
 713   where
 714    theta = (recip norms) NP.* (v1 `dot` v2)
 715    norms = (norm v1) NP.* (norm v2)
 716
 717
 718 -- | Retrieve the diagonal elements of the given matrix as a \"column
 719 --   vector,\" i.e. a m-by-1 matrix. We require the matrix to be
 720 --   square to avoid ambiguity in the return type which would ideally
 721 --   have dimension min(m,n) supposing an m-by-n matrix.
 722 --
 723 --   Examples:
 724 --
 725 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 726 --   >>> diagonal m
 727 --   ((1),(5),(9))
 728 --
 729 diagonal :: (Arity m) => Mat m m a -> Col m a
 730 diagonal matrix =
 731   construct lambda
 732   where
 733     lambda i _ = matrix !!! (i,i)
 734
 735
 736 -- | Given a square @matrix@, return a new matrix of the same size
 737 --   containing only the on-diagonal entries of @matrix@. The
 738 --   off-diagonal entries are set to zero.
 739 --
 740 --   Examples:
 741 --
 742 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 743 --   >>> diagonal_part m
 744 --   ((1,0,0),(0,5,0),(0,0,9))
 745 --
 746 diagonal_part :: (Arity m, Ring.C a)
 747          => Mat m m a
 748          -> Mat m m a
 749 diagonal_part matrix =
 750   construct lambda
 751   where
 752     lambda i j = if i == j then matrix !!! (i,j) else 0
 753
 754
 755 -- | Given a square @matrix@, return a new matrix of the same size
 756 --   containing only the on-diagonal and below-diagonal entries of
 757 --   @matrix@. The above-diagonal entries are set to zero.
 758 --
 759 --   Examples:
 760 --
 761 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 762 --   >>> lt_part m
 763 --   ((1,0,0),(4,5,0),(7,8,9))
 764 --
 765 lt_part :: (Arity m, Ring.C a)
 766         => Mat m m a
 767         -> Mat m m a
 768 lt_part matrix =
 769   construct lambda
 770   where
 771     lambda i j = if i >= j then matrix !!! (i,j) else 0
 772
 773
 774 -- | Given a square @matrix@, return a new matrix of the same size
 775 --   containing only the below-diagonal entries of @matrix@. The on-
 776 --   and above-diagonal entries are set to zero.
 777 --
 778 --   Examples:
 779 --
 780 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 781 --   >>> lt_part_strict m
 782 --   ((0,0,0),(4,0,0),(7,8,0))
 783 --
 784 lt_part_strict :: (Arity m, Ring.C a)
 785         => Mat m m a
 786         -> Mat m m a
 787 lt_part_strict matrix =
 788   construct lambda
 789   where
 790     lambda i j = if i > j then matrix !!! (i,j) else 0
 791
 792
 793 -- | Given a square @matrix@, return a new matrix of the same size
 794 --   containing only the on-diagonal and above-diagonal entries of
 795 --   @matrix@. The below-diagonal entries are set to zero.
 796 --
 797 --   Examples:
 798 --
 799 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 800 --   >>> ut_part m
 801 --   ((1,2,3),(0,5,6),(0,0,9))
 802 --
 803 ut_part :: (Arity m, Ring.C a)
 804         => Mat m m a
 805         -> Mat m m a
 806 ut_part = transpose . lt_part . transpose
 807
 808
 809 -- | Given a square @matrix@, return a new matrix of the same size
 810 --   containing only the above-diagonal entries of @matrix@. The on-
 811 --   and below-diagonal entries are set to zero.
 812 --
 813 --   Examples:
 814 --
 815 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 816 --   >>> ut_part_strict m
 817 --   ((0,2,3),(0,0,6),(0,0,0))
 818 --
 819 ut_part_strict :: (Arity m, Ring.C a)
 820         => Mat m m a
 821         -> Mat m m a
 822 ut_part_strict = transpose . lt_part_strict . transpose
 823
 824
 825 -- | Compute the trace of a square matrix, the sum of the elements
 826 --   which lie on its diagonal. We require the matrix to be
 827 --   square to avoid ambiguity in the return type which would ideally
 828 --   have dimension min(m,n) supposing an m-by-n matrix.
 829 --
 830 --   Examples:
 831 --
 832 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 833 --   >>> trace m
 834 --   15
 835 --
 836 trace :: (Arity m, Ring.C a) => Mat m m a -> a
 837 trace matrix =
 838   let (Mat rows) = diagonal matrix
 839   in
 840     element_sum $ V.map V.head rows
 841
 842
 843 -- | Zip together two matrices.
 844 --
 845 --   TODO: don't cheat with construct (map V.zips instead).
 846 --
 847 --   Examples:
 848 --
 849 --   >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
 850 --   >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
 851 --   >>> zip2 m1 m2
 852 --   (((1,1)),((1,2)),((1,3)))
 853 --
 854 --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
 855 --   >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
 856 --   >>> zip2 m1 m2
 857 --   (((1,1),(2,1)),((3,1),(4,1)))
 858 --
 859 zip2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a -> Mat m n (a,a)
 860 zip2 m1 m2 =
 861   construct lambda
 862   where
 863     lambda i j = (m1 !!! (i,j), m2 !!! (i,j))
 864
 865
 866 -- | Zip together three matrices.
 867 --
 868 --   TODO: don't cheat with construct (map V.zips instead).
 869 --
 870 --   Examples:
 871 --
 872 --   >>> let m1 = fromList [[1],[1],[1]] :: Col3 Int
 873 --   >>> let m2 = fromList [[1],[2],[3]] :: Col3 Int
 874 --   >>> let m3 = fromList [[4],[5],[6]] :: Col3 Int
 875 --   >>> zip2three m1 m2 m3
 876 --   (((1,1,4)),((1,2,5)),((1,3,6)))
 877 --
 878 --   >>> let m1 = fromList [[1,2],[3,4]] :: Mat2 Int
 879 --   >>> let m2 = fromList [[1,1],[1,1]] :: Mat2 Int
 880 --   >>> let m3 = fromList [[8,2],[6,3]] :: Mat2 Int
 881 --   >>> zip2three m1 m2 m3
 882 --   (((1,1,8),(2,1,2)),((3,1,6),(4,1,3)))
 883 --
 884 zip2three :: (Arity m, Arity n)
 885           => Mat m n a
 886           -> Mat m n a
 887           -> Mat m n a
 888           -> Mat m n (a,a,a)
 889 zip2three m1 m2 m3 =
 890   construct lambda
 891   where
 892     lambda i j = (m1 !!! (i,j), m2 !!! (i,j), m3 !!! (i,j))
 893
 894
 895 -- | Zip together two matrices using the supplied function.
 896 --
 897 --   Examples:
 898 --
 899 --   >>> let c1 = fromList [[1],[2],[3]] :: Col3 Integer
 900 --   >>> let c2 = fromList [[4],[5],[6]] :: Col3 Integer
 901 --   >>> zipwith2 (^) c1 c2
 902 --   ((1),(32),(729))
 903 --
 904 zipwith2 :: Arity m
 905            => (a -> a -> b)
 906            -> Col m a
 907            -> Col m a
 908            -> Col m b
 909 zipwith2 f c1 c2 =
 910   construct lambda
 911   where
 912     lambda i j = f (c1 !!! (i,j)) (c2 !!! (i,j))
 913
 914
 915 -- | Map a function over a matrix of any dimensions.
 916 --
 917 --   Examples:
 918 --
 919 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 920 --   >>> map2 (^2) m
 921 --   ((1,4),(9,16))
 922 --
 923 map2 :: (a -> b) -> Mat m n a -> Mat m n b
 924 map2 f (Mat rows) =
 925   Mat $ V.map g rows
 926   where
 927     g = V.map f
 928
 929
 930 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
 931 --   (of the row/column) to the accumulation function. The fold occurs
 932 --   from top-left to bottom-right.
 933 --
 934 --   Examples:
 935 --
 936 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 937 --   >>> ifoldl2 (\i j cur _ -> cur + i + j) 0 m
 938 --   18
 939 --
 940 ifoldl2 :: forall a b m n.
 941            (Int -> Int -> b -> a -> b)
 942         -> b
 943         -> Mat m n a
 944         -> b
 945 ifoldl2 f initial (Mat rows) =
 946   V.ifoldl row_function initial rows
 947   where
 948     -- | The order that we need this in (so that @g idx@ makes sense)
 949     --   is a little funny. So that we don't need to pass weird
 950     --   functions into ifoldl2, we swap the second and third
 951     --   arguments of @f@ calling the result @g@.
 952     g :: Int -> b -> Int -> a -> b
 953     g w x y = f w y x
 954
 955     row_function :: b -> Int -> Vec n a -> b
 956     row_function rowinit idx r = V.ifoldl (g idx) rowinit r
 957
 958
 959 -- | Fold over the entire matrix passing the coordinates @i@ and @j@
 960 --   (of the row/column) to the accumulation function. The fold occurs
 961 --   from bottom-right to top-left.
 962 --
 963 --   The order of the arguments in the supplied function are different
 964 --   from those in V.ifoldr; we keep them similar to ifoldl2.
 965 --
 966 --   Examples:
 967 --
 968 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
 969 --   >>> ifoldr2 (\i j cur _ -> cur + i + j) 0 m
 970 --   18
 971 --
 972 ifoldr2 :: forall a b m n.
 973            (Int -> Int -> b -> a -> b)
 974         -> b
 975         -> Mat m n a
 976         -> b
 977 ifoldr2 f initial (Mat rows) =
 978   V.ifoldr row_function initial rows
 979   where
 980     -- | Swap the order of arguments in @f@ so that it agrees with the
 981     --   @f@ passed to ifoldl2.
 982     g :: Int -> Int -> a -> b -> b
 983     g w x y z = f w x z y
 984
 985     row_function :: Int -> Vec n a -> b -> b
 986     row_function idx r rowinit = V.ifoldr (g idx) rowinit r
 987
 988
 989 -- | Map a function over a matrix of any dimensions, passing the
 990 --   coordinates @i@ and @j@ to the function @f@.
 991 --
 992 --   Examples:
 993 --
 994 --   >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
 995 --   >>> imap2 (\i j _ -> i+j) m
 996 --   ((0,1),(1,2))
 997 --
 998 imap2 :: (Int -> Int -> a -> b) -> Mat m n a -> Mat m n b
 999 imap2 f (Mat rows) =
1000   Mat $ V.imap g rows
1001   where
1002     g i = V.imap (f i)
1003
1004
1005 -- | Reverse the order of elements in a matrix.
1006 --
1007 --   Examples:
1008 --
1009 --   >>> let m1 = fromList [[1,2,3]] :: Row3 Int
1010 --   >>> reverse2 m1
1011 --   ((3,2,1))
1012 --
1013 --   >>> let m1 = vec3d (1,2,3 :: Int)
1014 --   >>> reverse2 m1
1015 --   ((3),(2),(1))
1016 --
1017 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1018 --   >>> reverse2 m
1019 --   ((9,8,7),(6,5,4),(3,2,1))
1020 --
1021 reverse2 :: (Arity m, Arity n) => Mat m n a -> Mat m n a
1022 reverse2 (Mat rows) = Mat $ V.reverse $ V.map V.reverse rows
1023
1024
1025 -- | Unsafely set the (i,j) element of the given matrix.
1026 --
1027 --   Examples:
1028 --
1029 --   >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
1030 --   >>> set_idx m (1,1) 17
1031 --   ((1,2,3),(4,17,6),(7,8,9))
1032 --
1033 set_idx :: forall m n a.
1034            (Arity m, Arity n)
1035         => Mat m n a
1036         -> (Int, Int)
1037         -> a
1038         -> Mat m n a
1039 set_idx matrix (i,j) newval =
1040   imap2 updater matrix
1041   where
1042     updater :: Int -> Int -> a -> a
1043     updater k l existing =
1044       if k == i && l == j
1045       then newval
1046       else existing
1047
1048
1049 -- | Compute the i,jth cofactor of the given @matrix@. This simply
1050 --   premultiplues the i,jth minor by (-1)^(i+j).
1051 cofactor :: (Arity m, Determined (Mat m m) a)
1052          => Mat (S m) (S m) a
1053          -> Int
1054          -> Int
1055          -> a
1056 cofactor matrix i j =
1057   (-1)^(toInteger i + toInteger j) NP.* (minor matrix i j)
1058
1059
1060 -- | Compute the inverse of a matrix using cofactor expansion
1061 --   (generalized Cramer's rule).
1062 --
1063 --   Examples:
1064 --
1065 --   >>> let m1 = fromList [[37,22],[17,54]] :: Mat2 Double
1066 --   >>> let e1 = [54/1624, -22/1624] :: [Double]
1067 --   >>> let e2 = [-17/1624, 37/1624] :: [Double]
1068 --   >>> let expected = fromList [e1, e2] :: Mat2 Double
1069 --   >>> let actual = inverse m1
1070 --   >>> frobenius_norm (actual - expected) < 1e-12
1071 --   True
1072 --
1073 inverse :: (Arity m,
1074             Determined (Mat (S m) (S m)) a,
1075             Determined (Mat m m) a,
1076             Field.C a)
1077         => Mat (S m) (S m) a
1078         -> Mat (S m) (S m) a
1079 inverse matrix =
1080   (1 / (determinant matrix)) *> (transpose $ construct lambda)
1081   where
1082     lambda i j = cofactor matrix i j
1083