{-# LANGUAGE RebindableSyntax #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{- |
Copyright    :   (c) Henning Thielemann 2009, Mikael Johansson 2006
Maintainer   :   numericprelude@henning-thielemann.de
Stability    :   provisional
Portability  :   requires multi-parameter type classes

Routines and abstractions for Matrices and
basic linear algebra over fields or rings.

We stick to simple Int indices.
Although advanced indices would be nice
e.g. for matrices with sub-matrices,
this is not easily implemented since arrays
do only support a lower and an upper bound
but no additional parameters.

ToDo:
 - Matrix inverse, determinant (see htam:Matrix)
-}

module MathObj.Matrix (
   T, Dimension,
   format,
   transpose,
   rows,
   columns,
   index,
   fromRows,
   fromColumns,
   fromList,
   dimension,
   numRows,
   numColumns,
   zipWith,
   zero,
   one,
   diagonal,
   scale,
   random,
   randomR,
   ) where

import qualified Algebra.Module   as Module
import qualified Algebra.Vector   as Vector
import qualified Algebra.Ring     as Ring
import qualified Algebra.Additive as Additive

import Algebra.Module((*>), )
import Algebra.Ring((*), fromInteger, scalarProduct, )
import Algebra.Additive((+), (-), subtract, )

import qualified System.Random as Rnd
import Data.Array (Array, array, listArray, accumArray, elems, bounds, (!), ixmap, range, )
import qualified Data.List as List

import Control.Monad (liftM2, )
import Control.Exception (assert, )

import Data.Function.HT (powerAssociative, )
import Data.Tuple.HT (swap, mapFst, )
import Data.List.HT (outerProduct, )
import Text.Show.HT (concatS, )

import NumericPrelude.Numeric (Int, )
import NumericPrelude.Base hiding (zipWith, )


{- $setup
>>> import qualified MathObj.Matrix as Matrix
>>> import qualified Algebra.Ring as Ring
>>> import qualified Algebra.Laws as Laws
>>> import Test.NumericPrelude.Utility ((/\))
>>> import qualified Test.QuickCheck as QC
>>> import NumericPrelude.Numeric as NP
>>> import NumericPrelude.Base as P
>>> import Prelude ()
>>>
>>> import Control.Monad (replicateM, join)
>>> import Control.Applicative (liftA2)
>>> import Data.Function.HT (nest)
>>>
>>> genDimension :: QC.Gen Int
>>> genDimension = QC.choose (0,20)
>>>
>>> genMatrixFor :: (QC.Arbitrary a) => Int -> Int -> QC.Gen (Matrix.T a)
>>> genMatrixFor m n =
>>>    fmap (Matrix.fromList m n) $ replicateM (m*n) QC.arbitrary
>>>
>>> genMatrix :: (QC.Arbitrary a) => QC.Gen (Matrix.T a)
>>> genMatrix = join $ liftA2 genMatrixFor genDimension genDimension
>>>
>>> genIntMatrix :: QC.Gen (Matrix.T Integer)
>>> genIntMatrix = genMatrix
>>>
>>> genFactorMatrix :: (QC.Arbitrary a) => Matrix.T a -> QC.Gen (Matrix.T a)
>>> genFactorMatrix a = genMatrixFor (Matrix.numColumns a) =<< genDimension
>>>
>>> genSameMatrix :: (QC.Arbitrary a) => Matrix.T a -> QC.Gen (Matrix.T a)
>>> genSameMatrix = uncurry genMatrixFor . Matrix.dimension
-}


{- |
A matrix is a twodimensional array, indexed by integers.
-}
newtype T a =
   Cons (Array (Dimension, Dimension) a)
      deriving (T a -> T a -> Bool
(T a -> T a -> Bool) -> (T a -> T a -> Bool) -> Eq (T a)
forall a. Eq a => T a -> T a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: T a -> T a -> Bool
$c/= :: forall a. Eq a => T a -> T a -> Bool
== :: T a -> T a -> Bool
$c== :: forall a. Eq a => T a -> T a -> Bool
Eq, Eq (T a)
Eq (T a)
-> (T a -> T a -> Ordering)
-> (T a -> T a -> Bool)
-> (T a -> T a -> Bool)
-> (T a -> T a -> Bool)
-> (T a -> T a -> Bool)
-> (T a -> T a -> T a)
-> (T a -> T a -> T a)
-> Ord (T a)
T a -> T a -> Bool
T a -> T a -> Ordering
T a -> T a -> T a
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
forall a. Ord a => Eq (T a)
forall a. Ord a => T a -> T a -> Bool
forall a. Ord a => T a -> T a -> Ordering
forall a. Ord a => T a -> T a -> T a
min :: T a -> T a -> T a
$cmin :: forall a. Ord a => T a -> T a -> T a
max :: T a -> T a -> T a
$cmax :: forall a. Ord a => T a -> T a -> T a
>= :: T a -> T a -> Bool
$c>= :: forall a. Ord a => T a -> T a -> Bool
> :: T a -> T a -> Bool
$c> :: forall a. Ord a => T a -> T a -> Bool
<= :: T a -> T a -> Bool
$c<= :: forall a. Ord a => T a -> T a -> Bool
< :: T a -> T a -> Bool
$c< :: forall a. Ord a => T a -> T a -> Bool
compare :: T a -> T a -> Ordering
$ccompare :: forall a. Ord a => T a -> T a -> Ordering
$cp1Ord :: forall a. Ord a => Eq (T a)
Ord, ReadPrec [T a]
ReadPrec (T a)
Int -> ReadS (T a)
ReadS [T a]
(Int -> ReadS (T a))
-> ReadS [T a] -> ReadPrec (T a) -> ReadPrec [T a] -> Read (T a)
forall a. Read a => ReadPrec [T a]
forall a. Read a => ReadPrec (T a)
forall a. Read a => Int -> ReadS (T a)
forall a. Read a => ReadS [T a]
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [T a]
$creadListPrec :: forall a. Read a => ReadPrec [T a]
readPrec :: ReadPrec (T a)
$creadPrec :: forall a. Read a => ReadPrec (T a)
readList :: ReadS [T a]
$creadList :: forall a. Read a => ReadS [T a]
readsPrec :: Int -> ReadS (T a)
$creadsPrec :: forall a. Read a => Int -> ReadS (T a)
Read)

type Dimension = Int

{- |
Transposition of matrices is just transposition in the sense of Data.List.

prop> genIntMatrix /\ \a -> Matrix.rows a == Matrix.columns (Matrix.transpose a)
prop> genIntMatrix /\ \a -> Matrix.columns a == Matrix.rows (Matrix.transpose a)
prop> genIntMatrix /\ \a -> genSameMatrix a /\ \b -> Laws.homomorphism Matrix.transpose (+) (+) a b
-}
transpose :: T a -> T a
transpose :: T a -> T a
transpose (Cons Array (Int, Int) a
m) =
   let ((Int, Int)
lower,(Int, Int)
upper) = Array (Int, Int) a -> ((Int, Int), (Int, Int))
forall i e. Array i e -> (i, i)
bounds Array (Int, Int) a
m
   in  Array (Int, Int) a -> T a
forall a. Array (Int, Int) a -> T a
Cons (((Int, Int), (Int, Int))
-> ((Int, Int) -> (Int, Int))
-> Array (Int, Int) a
-> Array (Int, Int) a
forall i j e.
(Ix i, Ix j) =>
(i, i) -> (i -> j) -> Array j e -> Array i e
ixmap ((Int, Int) -> (Int, Int)
forall a b. (a, b) -> (b, a)
swap (Int, Int)
lower, (Int, Int) -> (Int, Int)
forall a b. (a, b) -> (b, a)
swap (Int, Int)
upper) (Int, Int) -> (Int, Int)
forall a b. (a, b) -> (b, a)
swap Array (Int, Int) a
m)

rows :: T a -> [[a]]
rows :: T a -> [[a]]
rows mM :: T a
mM@(Cons Array (Int, Int) a
m) =
   let ((Int
lr,Int
lc), (Int
ur,Int
uc)) = Array (Int, Int) a -> ((Int, Int), (Int, Int))
forall i e. Array i e -> (i, i)
bounds Array (Int, Int) a
m
   in  (Int -> Int -> a) -> [Int] -> [Int] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [[c]]
outerProduct (T a -> Int -> Int -> a
forall a. T a -> Int -> Int -> a
index T a
mM) ((Int, Int) -> [Int]
forall a. Ix a => (a, a) -> [a]
range (Int
lr,Int
ur)) ((Int, Int) -> [Int]
forall a. Ix a => (a, a) -> [a]
range (Int
lc,Int
uc))

columns :: T a -> [[a]]
columns :: T a -> [[a]]
columns mM :: T a
mM@(Cons Array (Int, Int) a
m) =
   let ((Int
lr,Int
lc), (Int
ur,Int
uc)) = Array (Int, Int) a -> ((Int, Int), (Int, Int))
forall i e. Array i e -> (i, i)
bounds Array (Int, Int) a
m
   in  (Int -> Int -> a) -> [Int] -> [Int] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [[c]]
outerProduct ((Int -> Int -> a) -> Int -> Int -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip (T a -> Int -> Int -> a
forall a. T a -> Int -> Int -> a
index T a
mM)) ((Int, Int) -> [Int]
forall a. Ix a => (a, a) -> [a]
range (Int
lc,Int
uc)) ((Int, Int) -> [Int]
forall a. Ix a => (a, a) -> [a]
range (Int
lr,Int
ur))

index :: T a -> Dimension -> Dimension -> a
index :: T a -> Int -> Int -> a
index (Cons Array (Int, Int) a
m) Int
i Int
j = Array (Int, Int) a
m Array (Int, Int) a -> (Int, Int) -> a
forall i e. Ix i => Array i e -> i -> e
! (Int
i,Int
j)

{- |
prop> genIntMatrix /\ \a -> a == uncurry Matrix.fromRows (Matrix.dimension a) (Matrix.rows a)
-}
fromRows :: Dimension -> Dimension -> [[a]] -> T a
fromRows :: Int -> Int -> [[a]] -> T a
fromRows Int
m Int
n =
   Array (Int, Int) a -> T a
forall a. Array (Int, Int) a -> T a
Cons (Array (Int, Int) a -> T a)
-> ([[a]] -> Array (Int, Int) a) -> [[a]] -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ((Int, Int), (Int, Int)) -> [((Int, Int), a)] -> Array (Int, Int) a
forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int -> Int -> ((Int, Int), (Int, Int))
indexBounds Int
m Int
n) ([((Int, Int), a)] -> Array (Int, Int) a)
-> ([[a]] -> [((Int, Int), a)]) -> [[a]] -> Array (Int, Int) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   [[((Int, Int), a)]] -> [((Int, Int), a)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[((Int, Int), a)]] -> [((Int, Int), a)])
-> ([[a]] -> [[((Int, Int), a)]]) -> [[a]] -> [((Int, Int), a)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Int -> [(Int, a)] -> [((Int, Int), a)])
-> [Int] -> [[(Int, a)]] -> [[((Int, Int), a)]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
List.zipWith (\Int
r -> ((Int, a) -> ((Int, Int), a)) -> [(Int, a)] -> [((Int, Int), a)]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
c,a
x) -> ((Int
r,Int
c),a
x))) [Int]
allIndices ([[(Int, a)]] -> [[((Int, Int), a)]])
-> ([[a]] -> [[(Int, a)]]) -> [[a]] -> [[((Int, Int), a)]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ([a] -> [(Int, a)]) -> [[a]] -> [[(Int, a)]]
forall a b. (a -> b) -> [a] -> [b]
map ([Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
allIndices)

{- |
prop> genIntMatrix /\ \a -> a == uncurry Matrix.fromColumns (Matrix.dimension a) (Matrix.columns a)
-}
fromColumns :: Dimension -> Dimension -> [[a]] -> T a
fromColumns :: Int -> Int -> [[a]] -> T a
fromColumns Int
m Int
n =
   Array (Int, Int) a -> T a
forall a. Array (Int, Int) a -> T a
Cons (Array (Int, Int) a -> T a)
-> ([[a]] -> Array (Int, Int) a) -> [[a]] -> T a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ((Int, Int), (Int, Int)) -> [((Int, Int), a)] -> Array (Int, Int) a
forall i e. Ix i => (i, i) -> [(i, e)] -> Array i e
array (Int -> Int -> ((Int, Int), (Int, Int))
indexBounds Int
m Int
n) ([((Int, Int), a)] -> Array (Int, Int) a)
-> ([[a]] -> [((Int, Int), a)]) -> [[a]] -> Array (Int, Int) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   [[((Int, Int), a)]] -> [((Int, Int), a)]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[((Int, Int), a)]] -> [((Int, Int), a)])
-> ([[a]] -> [[((Int, Int), a)]]) -> [[a]] -> [((Int, Int), a)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Int -> [(Int, a)] -> [((Int, Int), a)])
-> [Int] -> [[(Int, a)]] -> [[((Int, Int), a)]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
List.zipWith (\Int
r -> ((Int, a) -> ((Int, Int), a)) -> [(Int, a)] -> [((Int, Int), a)]
forall a b. (a -> b) -> [a] -> [b]
map (\(Int
c,a
x) -> ((Int
c,Int
r),a
x))) [Int]
allIndices ([[(Int, a)]] -> [[((Int, Int), a)]])
-> ([[a]] -> [[(Int, a)]]) -> [[a]] -> [[((Int, Int), a)]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ([a] -> [(Int, a)]) -> [[a]] -> [[(Int, a)]]
forall a b. (a -> b) -> [a] -> [b]
map ([Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
allIndices)

fromList :: Dimension -> Dimension -> [a] -> T a
fromList :: Int -> Int -> [a] -> T a
fromList Int
m Int
n [a]
xs = Array (Int, Int) a -> T a
forall a. Array (Int, Int) a -> T a
Cons (((Int, Int), (Int, Int)) -> [a] -> Array (Int, Int) a
forall i e. Ix i => (i, i) -> [e] -> Array i e
listArray (Int -> Int -> ((Int, Int), (Int, Int))
indexBounds Int
m Int
n) [a]
xs)

appPrec :: Int
appPrec :: Int
appPrec = Int
10

instance (Show a) => Show (T a) where
   showsPrec :: Int -> T a -> ShowS
showsPrec Int
p T a
m =
      Bool -> ShowS -> ShowS
showParen (Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
appPrec)
         (String -> ShowS
showString String
"Matrix.fromRows " ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [[a]] -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
appPrec (T a -> [[a]]
forall a. T a -> [[a]]
rows T a
m))

format :: (Show a) => T a -> String
format :: T a -> String
format T a
m = T a -> ShowS
forall a. Show a => T a -> ShowS
formatS T a
m String
""

formatS :: (Show a) => T a -> ShowS
formatS :: T a -> ShowS
formatS =
   [ShowS] -> ShowS
concatS ([ShowS] -> ShowS) -> (T a -> [ShowS]) -> T a -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ([ShowS] -> ShowS) -> [[ShowS]] -> [ShowS]
forall a b. (a -> b) -> [a] -> [b]
map (\[ShowS]
r -> String -> ShowS
showString String
"(" ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [ShowS] -> ShowS
concatS [ShowS]
r ShowS -> ShowS -> ShowS
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> ShowS
showString String
")\n") ([[ShowS]] -> [ShowS]) -> (T a -> [[ShowS]]) -> T a -> [ShowS]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   ([a] -> [ShowS]) -> [[a]] -> [[ShowS]]
forall a b. (a -> b) -> [a] -> [b]
map (ShowS -> [ShowS] -> [ShowS]
forall a. a -> [a] -> [a]
List.intersperse (Char
' 'Char -> ShowS
forall a. a -> [a] -> [a]
:) ([ShowS] -> [ShowS]) -> ([a] -> [ShowS]) -> [a] -> [ShowS]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> ShowS) -> [a] -> [ShowS]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> a -> ShowS
forall a. Show a => Int -> a -> ShowS
showsPrec Int
11)) ([[a]] -> [[ShowS]]) -> (T a -> [[a]]) -> T a -> [[ShowS]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   T a -> [[a]]
forall a. T a -> [[a]]
rows

dimension :: T a -> (Dimension,Dimension)
dimension :: T a -> (Int, Int)
dimension (Cons Array (Int, Int) a
m) = ((Int, Int) -> (Int, Int) -> (Int, Int))
-> ((Int, Int), (Int, Int)) -> (Int, Int)
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry (Int, Int) -> (Int, Int) -> (Int, Int)
forall a. C a => a -> a -> a
subtract (Array (Int, Int) a -> ((Int, Int), (Int, Int))
forall i e. Array i e -> (i, i)
bounds Array (Int, Int) a
m) (Int, Int) -> (Int, Int) -> (Int, Int)
forall a. C a => a -> a -> a
+ (Int
1,Int
1)

numRows :: T a -> Dimension
numRows :: T a -> Int
numRows = (Int, Int) -> Int
forall a b. (a, b) -> a
fst ((Int, Int) -> Int) -> (T a -> (Int, Int)) -> T a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T a -> (Int, Int)
forall a. T a -> (Int, Int)
dimension

numColumns :: T a -> Dimension
numColumns :: T a -> Int
numColumns = (Int, Int) -> Int
forall a b. (a, b) -> b
snd ((Int, Int) -> Int) -> (T a -> (Int, Int)) -> T a -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. T a -> (Int, Int)
forall a. T a -> (Int, Int)
dimension

-- These implementations may benefit from a better exception than
-- just assertions to validate dimensionalities
{- |
prop> genIntMatrix /\ \a -> genSameMatrix a /\ \b -> Laws.commutative (+) a b
prop> genIntMatrix /\ \a -> genSameMatrix a /\ \b -> genSameMatrix b /\ \c -> Laws.associative (+) a b c
-}
instance (Additive.C a) => Additive.C (T a) where
   + :: T a -> T a -> T a
(+) = (a -> a -> a) -> T a -> T a -> T a
forall a b c. (a -> b -> c) -> T a -> T b -> T c
zipWith a -> a -> a
forall a. C a => a -> a -> a
(+)
   (-) = (a -> a -> a) -> T a -> T a -> T a
forall a b c. (a -> b -> c) -> T a -> T b -> T c
zipWith (-)
   zero :: T a
zero = Int -> Int -> T a
forall a. C a => Int -> Int -> T a
zero Int
1 Int
1

zipWith :: (a -> b -> c) -> T a -> T b -> T c
zipWith :: (a -> b -> c) -> T a -> T b -> T c
zipWith a -> b -> c
op mM :: T a
mM@(Cons Array (Int, Int) a
m) nM :: T b
nM@(Cons Array (Int, Int) b
n) =
   let d :: (Int, Int)
d = T a -> (Int, Int)
forall a. T a -> (Int, Int)
dimension T a
mM
       em :: [a]
em = Array (Int, Int) a -> [a]
forall i e. Array i e -> [e]
elems Array (Int, Int) a
m
       en :: [b]
en = Array (Int, Int) b -> [b]
forall i e. Array i e -> [e]
elems Array (Int, Int) b
n
   in  Bool -> T c -> T c
forall a. (?callStack::CallStack) => Bool -> a -> a
assert ((Int, Int)
d (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
== T b -> (Int, Int)
forall a. T a -> (Int, Int)
dimension T b
nM) (T c -> T c) -> T c -> T c
forall a b. (a -> b) -> a -> b
$
         (Int -> Int -> [c] -> T c) -> (Int, Int) -> [c] -> T c
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Int -> Int -> [c] -> T c
forall a. Int -> Int -> [a] -> T a
fromList (Int, Int)
d ((a -> b -> c) -> [a] -> [b] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
List.zipWith a -> b -> c
op [a]
em [b]
en)

{- |
prop> genIntMatrix /\ \a -> Laws.identity (+) (uncurry Matrix.zero $ Matrix.dimension a) a
-}
zero :: (Additive.C a) => Dimension -> Dimension -> T a
zero :: Int -> Int -> T a
zero Int
m Int
n =
   Int -> Int -> [a] -> T a
forall a. Int -> Int -> [a] -> T a
fromList Int
m Int
n ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$
   a -> [a]
forall a. a -> [a]
List.repeat a
forall a. C a => a
Additive.zero
--    List.replicate (fromInteger (m*n)) zero

one :: (Ring.C a) => Dimension -> T a
one :: Int -> T a
one Int
n =
   Array (Int, Int) a -> T a
forall a. Array (Int, Int) a -> T a
Cons (Array (Int, Int) a -> T a) -> Array (Int, Int) a -> T a
forall a b. (a -> b) -> a -> b
$
   (a -> a -> a)
-> a
-> ((Int, Int), (Int, Int))
-> [((Int, Int), a)]
-> Array (Int, Int) a
forall i e a.
Ix i =>
(e -> a -> e) -> e -> (i, i) -> [(i, a)] -> Array i e
accumArray ((a -> a -> a) -> a -> a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> a -> a
forall a b. a -> b -> a
const) a
forall a. C a => a
Additive.zero
      (Int -> Int -> ((Int, Int), (Int, Int))
indexBounds Int
n Int
n)
      ((Int -> ((Int, Int), a)) -> [Int] -> [((Int, Int), a)]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
i -> ((Int
i,Int
i), a
forall a. C a => a
Ring.one)) (Int -> [Int]
indexRange Int
n))

{- |
prop> genDimension /\ \n -> Matrix.one n == Matrix.diagonal (replicate n Ring.one :: [Integer])
-}
diagonal :: (Additive.C a) => [a] -> T a
diagonal :: [a] -> T a
diagonal [a]
xs =
   let n :: Int
n = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
List.length [a]
xs
   in  Array (Int, Int) a -> T a
forall a. Array (Int, Int) a -> T a
Cons (Array (Int, Int) a -> T a) -> Array (Int, Int) a -> T a
forall a b. (a -> b) -> a -> b
$
       (a -> a -> a)
-> a
-> ((Int, Int), (Int, Int))
-> [((Int, Int), a)]
-> Array (Int, Int) a
forall i e a.
Ix i =>
(e -> a -> e) -> e -> (i, i) -> [(i, a)] -> Array i e
accumArray ((a -> a -> a) -> a -> a -> a
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> a -> a
forall a b. a -> b -> a
const) a
forall a. C a => a
Additive.zero
          (Int -> Int -> ((Int, Int), (Int, Int))
indexBounds Int
n Int
n)
          ([(Int, Int)] -> [a] -> [((Int, Int), a)]
forall a b. [a] -> [b] -> [(a, b)]
zip ((Int -> (Int, Int)) -> [Int] -> [(Int, Int)]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
i -> (Int
i,Int
i)) [Int]
allIndices) [a]
xs)

scale :: (Ring.C a) => a -> T a -> T a
scale :: a -> T a -> T a
scale a
s = a -> T a -> T a
forall (v :: * -> *) a. (Functor v, C a) => a -> v a -> v a
Vector.functorScale a
s

{- |
prop> genIntMatrix /\ \a -> Laws.leftIdentity  (*) (Matrix.one (Matrix.numRows a)) a
prop> genIntMatrix /\ \a -> Laws.rightIdentity (*) (Matrix.one (Matrix.numColumns a)) a
prop> genIntMatrix /\ \a -> genFactorMatrix a /\ \b -> Laws.homomorphism Matrix.transpose (*) (flip (*)) a b
prop> genIntMatrix /\ \a -> genFactorMatrix a /\ \b -> genFactorMatrix b /\ \c -> Laws.associative (*) a b c
prop> genIntMatrix /\ \b -> genSameMatrix b /\ \c -> genFactorMatrix b /\ \a -> Laws.leftDistributive (*) (+) a b c
prop> genIntMatrix /\ \a -> genFactorMatrix a /\ \b -> genSameMatrix b /\ \c -> Laws.rightDistributive (*) (+) a b c
prop> QC.choose (0,10) /\ \k -> genDimension /\ \n -> genMatrixFor n n /\ \a -> a^k == nest (fromInteger k) ((a::Matrix.T Integer)*) (Matrix.one n)
-}
instance (Ring.C a) => Ring.C (T a) where
   T a
mM * :: T a -> T a -> T a
* T a
nM =
      Bool -> T a -> T a
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (T a -> Int
forall a. T a -> Int
numColumns T a
mM Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== T a -> Int
forall a. T a -> Int
numRows T a
nM) (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$
      Int -> Int -> [a] -> T a
forall a. Int -> Int -> [a] -> T a
fromList (T a -> Int
forall a. T a -> Int
numRows T a
mM) (T a -> Int
forall a. T a -> Int
numColumns T a
nM) ([a] -> T a) -> [a] -> T a
forall a b. (a -> b) -> a -> b
$
      ([a] -> [a] -> a) -> [[a]] -> [[a]] -> [a]
forall (m :: * -> *) a1 a2 r.
Monad m =>
(a1 -> a2 -> r) -> m a1 -> m a2 -> m r
liftM2 [a] -> [a] -> a
forall a. C a => [a] -> [a] -> a
scalarProduct (T a -> [[a]]
forall a. T a -> [[a]]
rows T a
mM) (T a -> [[a]]
forall a. T a -> [[a]]
columns T a
nM)
   fromInteger :: Integer -> T a
fromInteger Integer
n = Int -> Int -> [a] -> T a
forall a. Int -> Int -> [a] -> T a
fromList Int
1 Int
1 [Integer -> a
forall a. C a => Integer -> a
fromInteger Integer
n]
   T a
mM ^ :: T a -> Integer -> T a
^ Integer
n =
      Bool -> T a -> T a
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (T a -> Int
forall a. T a -> Int
numColumns T a
mM Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== T a -> Int
forall a. T a -> Int
numRows T a
mM) (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$
      Bool -> T a -> T a
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
forall a. C a => a
Additive.zero) (T a -> T a) -> T a -> T a
forall a b. (a -> b) -> a -> b
$
      (T a -> T a -> T a) -> T a -> T a -> Integer -> T a
forall a. (a -> a -> a) -> a -> a -> Integer -> a
powerAssociative T a -> T a -> T a
forall a. C a => a -> a -> a
(*) (Int -> T a
forall a. C a => Int -> T a
one (T a -> Int
forall a. T a -> Int
numColumns T a
mM)) T a
mM Integer
n

instance Functor T where
   fmap :: (a -> b) -> T a -> T b
fmap a -> b
f (Cons Array (Int, Int) a
m) = Array (Int, Int) b -> T b
forall a. Array (Int, Int) a -> T a
Cons ((a -> b) -> Array (Int, Int) a -> Array (Int, Int) b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> b
f Array (Int, Int) a
m)

instance Vector.C T where
   zero :: T a
zero  = T a
forall a. C a => a
Additive.zero
   <+> :: T a -> T a -> T a
(<+>) = T a -> T a -> T a
forall a. C a => a -> a -> a
(+)
   *> :: a -> T a -> T a
(*>)  = a -> T a -> T a
forall a. C a => a -> T a -> T a
scale

instance Module.C a b => Module.C a (T b) where
   a
x *> :: a -> T b -> T b
*> T b
m = (b -> b) -> T b -> T b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (a
xa -> b -> b
forall a v. C a v => a -> v -> v
*>) T b
m


random :: (Rnd.RandomGen g, Rnd.Random a) =>
   Dimension -> Dimension -> g -> (T a, g)
random :: Int -> Int -> g -> (T a, g)
random =
   (g -> (a, g)) -> Int -> Int -> g -> (T a, g)
forall g a.
(RandomGen g, Random a) =>
(g -> (a, g)) -> Int -> Int -> g -> (T a, g)
randomAux g -> (a, g)
forall a g. (Random a, RandomGen g) => g -> (a, g)
Rnd.random

randomR :: (Rnd.RandomGen g, Rnd.Random a) =>
   Dimension -> Dimension -> (a,a) -> g -> (T a, g)
randomR :: Int -> Int -> (a, a) -> g -> (T a, g)
randomR Int
m Int
n (a, a)
rng =
   (g -> (a, g)) -> Int -> Int -> g -> (T a, g)
forall g a.
(RandomGen g, Random a) =>
(g -> (a, g)) -> Int -> Int -> g -> (T a, g)
randomAux ((a, a) -> g -> (a, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
Rnd.randomR (a, a)
rng) Int
m Int
n

{-
could be made nicer with the State monad,
but I like to keep dependencies minimal
-}
randomAux :: (Rnd.RandomGen g, Rnd.Random a) =>
   (g -> (a,g)) -> Dimension -> Dimension -> g -> (T a, g)
randomAux :: (g -> (a, g)) -> Int -> Int -> g -> (T a, g)
randomAux g -> (a, g)
rnd Int
m Int
n g
g0 =
   ([a] -> T a) -> ([a], g) -> (T a, g)
forall a c b. (a -> c) -> (a, b) -> (c, b)
mapFst (Int -> Int -> [a] -> T a
forall a. Int -> Int -> [a] -> T a
fromList Int
m Int
n) (([a], g) -> (T a, g)) -> ([a], g) -> (T a, g)
forall a b. (a -> b) -> a -> b
$ (g, [a]) -> ([a], g)
forall a b. (a, b) -> (b, a)
swap ((g, [a]) -> ([a], g)) -> (g, [a]) -> ([a], g)
forall a b. (a -> b) -> a -> b
$
   (g -> Int -> (g, a)) -> g -> [Int] -> (g, [a])
forall (t :: * -> *) a b c.
Traversable t =>
(a -> b -> (a, c)) -> a -> t b -> (a, t c)
List.mapAccumL (\g
g Int
_i -> (a, g) -> (g, a)
forall a b. (a, b) -> (b, a)
swap ((a, g) -> (g, a)) -> (a, g) -> (g, a)
forall a b. (a -> b) -> a -> b
$ g -> (a, g)
rnd g
g) g
g0 (Int -> [Int]
indexRange (Int
mInt -> Int -> Int
forall a. C a => a -> a -> a
*Int
n))

{-
What more do we need from our matrix type? We have addition,
subtraction and multiplication, and thus composition of generic
free-module-maps. We're going to want to solve linear equations with
or without fields underneath, so we're going to want an implementation
of the Gaussian algorithm as well as most probably Smith normal
form. Determinants are cool, and these are to be calculated either
with the Gaussian algorithm or some other goodish method.
-}

{-
{- |
 We'll want generic linear equation solving, returning one solution,
any solution really, or nothing. Basically, this is asking for the
preimage of a given vector over the given map, so

a_11 x_1 + .. + a_1n x_n = y_1
 ...
a_m1 x_1 + .. + a_mn a_n = y_m

has really x_1,...,x_n as a preimage of the vector y_1,..,y_m under
the map (a_ij), since obviously y_1,..,y_m = (a_ij) x_1,..,x_n

So, generic linear equation solving boils down to the function
-}
preimage :: (Ring.C a) => T a -> T a -> Maybe (T a)
preimage a y = assert
        (numRows a == numRows y &&     -- they match
         numColumns y == 1)               -- and y is a column vector
                Nothing
-}

{-
Cf. /usr/lib/hugs/demos/Matrix.hs
-}


-- these functions control whether we use 0 or 1 based indices

indexRange :: Dimension -> [Dimension]
indexRange :: Int -> [Int]
indexRange Int
n = [Int
0..(Int
nInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1)]

indexBounds ::
   Dimension -> Dimension ->
   ((Dimension,Dimension), (Dimension,Dimension))
indexBounds :: Int -> Int -> ((Int, Int), (Int, Int))
indexBounds Int
m Int
n =
   ((Int
0,Int
0), (Int
mInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1,Int
nInt -> Int -> Int
forall a. C a => a -> a -> a
-Int
1))

allIndices :: [Dimension]
allIndices :: [Int]
allIndices = [Int
0..]