module AERN2.Linear.Matrix.Type where

import qualified Prelude as P
import MixedTypesNumPrelude
import qualified Numeric.CollectErrors as CN
import AERN2.Linear.Vector.Type (Vector, (!))
import qualified AERN2.Linear.Vector.Type as V
import qualified Data.List as List

import AERN2.MP.Ball

data (Matrix a) = 
    Matrix
    {
        Matrix a -> Integer
width     :: Integer,
        Matrix a -> Vector a
entries   :: Vector a
    } deriving (Int -> Matrix a -> ShowS
[Matrix a] -> ShowS
Matrix a -> String
(Int -> Matrix a -> ShowS)
-> (Matrix a -> String) -> ([Matrix a] -> ShowS) -> Show (Matrix a)
forall a. Show a => Int -> Matrix a -> ShowS
forall a. Show a => [Matrix a] -> ShowS
forall a. Show a => Matrix a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Matrix a] -> ShowS
$cshowList :: forall a. Show a => [Matrix a] -> ShowS
show :: Matrix a -> String
$cshow :: forall a. Show a => Matrix a -> String
showsPrec :: Int -> Matrix a -> ShowS
$cshowsPrec :: forall a. Show a => Int -> Matrix a -> ShowS
Show)

height :: Matrix a -> Integer
height :: Matrix a -> Integer
height (Matrix Integer
w Vector a
e) = 
    (Vector a -> Integer
forall a. Vector a -> Integer
V.length Vector a
e) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`P.div` Integer
w

get :: Matrix a -> Integer -> Integer -> a
get :: Matrix a -> Integer -> Integer -> a
get Matrix a
m Integer
i Integer
j =
    Matrix a -> Vector a
forall a. Matrix a -> Vector a
entries Matrix a
m Vector a -> Integer -> a
forall a. Vector a -> Integer -> a
! (Integer
i Integer -> Integer -> MulType Integer Integer
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* (Matrix a -> Integer
forall a. Matrix a -> Integer
width Matrix a
m) Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
j)

identity :: (HasIntegers a) => Integer -> Integer -> Matrix a
identity :: Integer -> Integer -> Matrix a
identity Integer
m Integer
n = 
    Integer -> Integer -> a -> Matrix a
forall a. HasIntegers a => Integer -> Integer -> a -> Matrix a
diag Integer
m Integer
n (Integer -> a
forall t1 t2. ConvertibleExactly t1 t2 => t1 -> t2
convertExactly Integer
1)

diag :: (HasIntegers a) => Integer -> Integer -> a -> Matrix a
diag :: Integer -> Integer -> a -> Matrix a
diag Integer
m Integer
n a
x = 
    Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
forall a.
Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
create Integer
m Integer
n (\Integer
i Integer
j -> if Integer
i Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
j then a
x else (Integer -> a
forall t1 t2. ConvertibleExactly t1 t2 => t1 -> t2
convertExactly Integer
0))

rows :: Matrix a -> [Vector a]
rows :: Matrix a -> [Vector a]
rows m :: Matrix a
m@(Matrix Integer
w Vector a
e) = 
    [Integer -> Integer -> Vector a -> Vector a
forall a. Integer -> Integer -> Vector a -> Vector a
V.slice (Integer
iInteger -> Integer -> MulType Integer Integer
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
*Integer
w) Integer
w Vector a
e| Integer
i <- [Integer
0 .. Matrix a -> Integer
forall a. Matrix a -> Integer
height Matrix a
m Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1]]

columns :: Matrix a -> Vector (Vector a)
columns :: Matrix a -> Vector (Vector a)
columns Matrix a
m = 
    (Integer -> Vector a) -> Vector Integer -> Vector (Vector a)
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\Integer
j -> (Integer -> a) -> Vector Integer -> Vector a
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\Integer
i -> Matrix a -> Integer -> Integer -> a
forall a. Matrix a -> Integer -> Integer -> a
get Matrix a
m Integer
i Integer
j) (Vector Integer -> Vector a) -> Vector Integer -> Vector a
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Vector Integer
forall a. Enum a => a -> a -> Vector a
V.enumFromTo Integer
0 (Matrix a -> Integer
forall a. Matrix a -> Integer
height Matrix a
m Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1)) (Vector Integer -> Vector (Vector a))
-> Vector Integer -> Vector (Vector a)
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Vector Integer
forall a. Enum a => a -> a -> Vector a
V.enumFromTo Integer
0 (Matrix a -> Integer
forall a. Matrix a -> Integer
width Matrix a
m Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1)

create :: Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
create :: Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
create Integer
m Integer
n Integer -> Integer -> a
f =
    Integer -> Vector a -> Matrix a
forall a. Integer -> Vector a -> Matrix a
Matrix Integer
n (Vector a -> Matrix a) -> Vector a -> Matrix a
forall a b. (a -> b) -> a -> b
$ (Integer -> a) -> Vector Integer -> Vector a
forall a b. (a -> b) -> Vector a -> Vector b
V.map (\Integer
x -> Integer -> Integer -> a
f (Integer -> Integer
i Integer
x) (Integer -> ModType Integer Integer
j Integer
x)) (Vector Integer -> Vector a) -> Vector Integer -> Vector a
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Vector Integer
forall a. Enum a => a -> a -> Vector a
V.enumFromTo Integer
0 (Integer
mInteger -> Integer -> MulType Integer Integer
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
*Integer
n Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1)
    where
    j :: Integer -> ModType Integer Integer
j Integer
x = Integer
x Integer -> Integer -> ModType Integer Integer
forall t1 t2. CanDivIMod t1 t2 => t1 -> t2 -> ModType t1 t2
`mod` Integer
n
    i :: Integer -> Integer
i Integer
x = (Integer
x Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer -> ModType Integer Integer
j Integer
x) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`P.div` Integer
n

imap :: (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a
imap :: (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a
imap Integer -> Integer -> a -> a
f (Matrix Integer
w Vector a
ents) =
    Integer -> Vector a -> Matrix a
forall a. Integer -> Vector a -> Matrix a
Matrix Integer
w ((Integer -> a -> a) -> Vector a -> Vector a
forall a b. (Integer -> a -> b) -> Vector a -> Vector b
V.imap Integer -> a -> a
g Vector a
ents)
    where
    j :: Integer -> ModType Integer Integer
j Integer
x = Integer
x Integer -> Integer -> ModType Integer Integer
forall t1 t2. CanDivIMod t1 t2 => t1 -> t2 -> ModType t1 t2
`mod` Integer
w
    i :: Integer -> Integer
i Integer
x = (Integer
x Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer -> ModType Integer Integer
j Integer
x) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`P.div` Integer
w
    g :: Integer -> a -> a
g Integer
k a
x = Integer -> Integer -> a -> a
f (Integer -> Integer
i Integer
k) (Integer -> ModType Integer Integer
j Integer
k) a
x

instance CanIntersectAsymmetric (Matrix (CN MPBall)) (Matrix (CN MPBall)) where
    type IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall)) = Matrix (CN MPBall)
    intersect :: Matrix (CN MPBall)
-> Matrix (CN MPBall)
-> IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall))
intersect (Matrix Integer
w0 Vector (CN MPBall)
v0) (Matrix Integer
_w1 Vector (CN MPBall)
v1) =
        Integer -> Vector (CN MPBall) -> Matrix (CN MPBall)
forall a. Integer -> Vector a -> Matrix a
Matrix Integer
w0 (Vector (CN MPBall) -> Matrix (CN MPBall))
-> Vector (CN MPBall) -> Matrix (CN MPBall)
forall a b. (a -> b) -> a -> b
$ (CN MPBall -> CN MPBall -> CN MPBall)
-> Vector (CN MPBall) -> Vector (CN MPBall) -> Vector (CN MPBall)
forall a b c. (a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith CN MPBall -> CN MPBall -> CN MPBall
forall e1 e2.
CanIntersectAsymmetric e1 e2 =>
e1 -> e2 -> IntersectionType e1 e2
intersect Vector (CN MPBall)
v0 Vector (CN MPBall)
v1

inftyNorm :: (CanAddSameType a, CanSubSameType a, CanAbsSameType a, HasIntegers a, CanMinMaxSameType a) => Matrix a -> a
inftyNorm :: Matrix a -> a
inftyNorm (Matrix a
m :: Matrix a) = 
    -- TODO: could be optimised.
    (a -> a -> a) -> a -> [a] -> a
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
List.foldl' a -> a -> a
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max (Integer -> a
forall t1 t2. ConvertibleExactly t1 t2 => t1 -> t2
convertExactly Integer
0 :: a)
    [
        (a -> a -> a) -> a -> Vector a -> a
forall b a. (b -> a -> b) -> b -> Vector a -> b
V.foldl' a -> a -> a
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
(+) (Integer -> a
forall t1 t2. ConvertibleExactly t1 t2 => t1 -> t2
convertExactly Integer
0 :: a) (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a) -> Vector a -> Vector a
forall a b. (a -> b) -> Vector a -> Vector b
V.map a -> a
forall t. CanAbs t => t -> AbsType t
abs Vector a
r 
        |
        Vector a
r <- Matrix a -> [Vector a]
forall a. Matrix a -> [Vector a]
rows Matrix a
m
    ]

instance Functor Matrix where
    fmap :: (a -> b) -> Matrix a -> Matrix b
fmap a -> b
h Matrix a
m = 
        Integer -> Vector b -> Matrix b
forall a. Integer -> Vector a -> Matrix a
Matrix (Matrix a -> Integer
forall a. Matrix a -> Integer
width Matrix a
m) ((a -> b) -> Vector a -> Vector b
forall a b. (a -> b) -> Vector a -> Vector b
V.map a -> b
h (Matrix a -> Vector a
forall a. Matrix a -> Vector a
entries Matrix a
m))

instance 
    (CanAddSameType a, CanMulSameType a, HasIntegers a) =>
    CanMulAsymmetric (Matrix a) (Matrix a)
    where
    type MulType (Matrix a) (Matrix a) = Matrix a
    mul :: Matrix a -> Matrix a -> MulType (Matrix a) (Matrix a)
mul Matrix a
m0 Matrix a
m1 = 
        Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
forall a.
Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
create (Matrix a -> Integer
forall a. Matrix a -> Integer
height Matrix a
m0) (Matrix a -> Integer
forall a. Matrix a -> Integer
width Matrix a
m1) (Integer -> a -> Integer -> Integer -> a
aux Integer
0 (Integer -> a
forall t1 t2. ConvertibleExactly t1 t2 => t1 -> t2
convertExactly Integer
0))
        where
        aux :: Integer -> a -> Integer -> Integer -> a
aux Integer
k a
sm Integer
i Integer
j = 
            if Integer
k Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Matrix a -> Integer
forall a. Matrix a -> Integer
width Matrix a
m0 then 
                a
sm
            else 
                Integer -> a -> Integer -> Integer -> a
aux (Integer
k Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
1) (a
sm a -> a -> AddType a a
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ (Matrix a -> Integer -> Integer -> a
forall a. Matrix a -> Integer -> Integer -> a
get Matrix a
m0 Integer
i Integer
k) a -> a -> MulType a a
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* (Matrix a -> Integer -> Integer -> a
forall a. Matrix a -> Integer -> Integer -> a
get Matrix a
m1 Integer
k Integer
j)) Integer
i Integer
j


instance 
    (CanAddSameType a) =>
    CanAddAsymmetric (Matrix a) (Matrix a)
    where
    type AddType (Matrix a) (Matrix a) = Matrix a
    add :: Matrix a -> Matrix a -> AddType (Matrix a) (Matrix a)
add (Matrix Integer
w Vector a
e) (Matrix Integer
_ Vector a
e') =
        Integer -> Vector a -> Matrix a
forall a. Integer -> Vector a -> Matrix a
Matrix Integer
w (Vector a
e Vector a -> Vector a -> AddType (Vector a) (Vector a)
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Vector a
e')

instance 
    (CanSubSameType a) =>
    CanSub (Matrix a) (Matrix a)
    where
    type SubType (Matrix a) (Matrix a) = Matrix a
    sub :: Matrix a -> Matrix a -> SubType (Matrix a) (Matrix a)
sub (Matrix Integer
w Vector a
e) (Matrix Integer
_ Vector a
e') =
        Integer -> Vector a -> Matrix a
forall a. Integer -> Vector a -> Matrix a
Matrix Integer
w (Vector a
e Vector a -> Vector a -> SubType (Vector a) (Vector a)
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Vector a
e')
    

instance 
    (CanAddSameType a, CanMulSameType a, HasIntegers a) =>
    CanMulAsymmetric (Matrix a) (Vector a)
    where
    type MulType (Matrix a) (Vector a) = Vector a
    mul :: Matrix a -> Vector a -> MulType (Matrix a) (Vector a)
mul m :: Matrix a
m@(Matrix Integer
_w Vector a
_e) Vector a
v =
        [a] -> Vector a
forall a. [a] -> Vector a
V.fromList [Vector a
r Vector a -> Vector a -> MulType (Vector a) (Vector a)
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Vector a
v| Vector a
r <- Matrix a -> [Vector a]
forall a. Matrix a -> [Vector a]
rows Matrix a
m]

instance 
    (HasAccuracy a, HasPrecision a) => HasAccuracy (Matrix a)
    where
    getAccuracy :: Matrix a -> Accuracy
getAccuracy Matrix a
m = 
        (Accuracy -> Accuracy -> Accuracy)
-> Accuracy -> Vector Accuracy -> Accuracy
forall b a. (b -> a -> b) -> b -> Vector a -> b
V.foldl' Accuracy -> Accuracy -> Accuracy
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max Accuracy
NoInformation (Vector Accuracy -> Accuracy) -> Vector Accuracy -> Accuracy
forall a b. (a -> b) -> a -> b
$ (a -> Accuracy) -> Vector a -> Vector Accuracy
forall a b. (a -> b) -> Vector a -> Vector b
V.map a -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy (Matrix a -> Vector a
forall a. Matrix a -> Vector a
entries Matrix a
m)

instance 
    (HasPrecision a) => HasPrecision (Matrix a)
    where
    getPrecision :: Matrix a -> Precision
getPrecision Matrix a
m = 
        (Precision -> Precision -> Precision)
-> Precision -> Vector Precision -> Precision
forall b a. (b -> a -> b) -> b -> Vector a -> b
V.foldl' Precision -> Precision -> Precision
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max (Integer -> Precision
prec Integer
2) (Vector Precision -> Precision) -> Vector Precision -> Precision
forall a b. (a -> b) -> a -> b
$ (a -> Precision) -> Vector a -> Vector Precision
forall a b. (a -> b) -> Vector a -> Vector b
V.map a -> Precision
forall t. HasPrecision t => t -> Precision
getPrecision (Matrix a -> Vector a
forall a. Matrix a -> Vector a
entries Matrix a
m)

instance 
    (CN.CanTestErrorsPresent a) => CN.CanTestErrorsPresent (Matrix a)
    where
    hasError :: Matrix a -> Bool
hasError Matrix a
m = (Bool -> Bool -> Bool) -> Bool -> Vector Bool -> Bool
forall b a. (b -> a -> b) -> b -> Vector a -> b
V.foldl' Bool -> Bool -> Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
(||) Bool
False (Vector Bool -> Bool) -> Vector Bool -> Bool
forall a b. (a -> b) -> a -> b
$ (a -> Bool) -> Vector a -> Vector Bool
forall a b. (a -> b) -> Vector a -> Vector b
V.map (a -> Bool
forall es. CanTestErrorsPresent es => es -> Bool
CN.hasError) (Matrix a -> Vector a
forall a. Matrix a -> Vector a
entries Matrix a
m)