-- SPDX-License-Identifier: MPL-2.0
module HMath.Matrix.Internal (
    entries,
    build,
    get,
    mapEntries,
    finiteMatrix,
    finiteValues,
    stableNorm,
    dotV,
    replace,
    swap,
    epsilon,
    defaultTolerance,
) where

import HMath.Error
import HMath.Matrix

entries :: Matrix a -> [[a]]
entries :: forall a. Matrix a -> [[a]]
entries Matrix a
a = let (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a in [Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
c (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
c) (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a)) | Int
i <- [Int
0 .. Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]

build :: Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build :: forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
r Int
c Int -> Int -> a
f =
    Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked
        Int
r
        Int
c
        (if Int
r Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
c Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then [] else [Int -> Int -> a
f Int
i Int
j | Int
i <- [Int
0 .. Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
j <- [Int
0 .. Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]])

get :: Matrix a -> Int -> Int -> a
get :: forall a. Matrix a -> Int -> Int -> a
get Matrix a
a Int
i Int
j = Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j)

mapEntries :: (a -> b) -> Matrix a -> Matrix b
mapEntries :: forall a b. (a -> b) -> Matrix a -> Matrix b
mapEntries a -> b
f Matrix a
a = case Int -> Int -> [b] -> Either MathError (Matrix b)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
r Int
c ((a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map a -> b
f (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a)) of
    Right Matrix b
b -> Matrix b
b
    Left MathError
_ -> [Char] -> Matrix b
forall a. HasCallStack => [Char] -> a
error [Char]
"mapEntries: internal shape invariant"
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a

finiteMatrix :: (RealFloat a) => Matrix a -> Either MathError ()
finiteMatrix :: forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix a
a = case [Int
k | (Int
k, a
x) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a), a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x] of
    Int
k : [Int]
_ -> MathError -> Either MathError ()
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonFiniteMatrixEntry (Int
k Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a), Int
k Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a)))
    [] -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()

finiteValues :: (RealFloat a) => [a] -> Either MathError ()
finiteValues :: forall a. RealFloat a => [a] -> Either MathError ()
finiteValues [a]
xs = case [Int
i | (Int
i, a
x) <- [Int] -> [a] -> [(Int, a)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
0 ..] [a]
xs, a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
x Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
x] of
    Int
i : [Int]
_ -> MathError -> Either MathError ()
forall a b. a -> Either a b
Left (Int -> MathError
NonFiniteVectorEntry Int
i)
    [] -> () -> Either MathError ()
forall a b. b -> Either a b
Right ()

stableNorm :: (RealFloat a) => [a] -> a
stableNorm :: forall a. RealFloat a => [a] -> a
stableNorm [a]
xs
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any a -> Bool
forall a. RealFloat a => a -> Bool
isNaN [a]
xs = a
0 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0
    | (a -> Bool) -> [a] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite [a]
xs = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
0
    | a
s a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a
0
    | Bool
otherwise = a
s a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sqrt ([a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [(a
x a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
s) a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) | a
x <- [a]
xs])
  where
    s :: a
s = [a] -> a
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (a
0 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
forall a. Num a => a -> a
abs [a]
xs)

dotV :: (Num a) => [a] -> [a] -> a
dotV :: forall a. Num a => [a] -> [a] -> a
dotV [a]
a [a]
b = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
a [a]
b)

replace :: Int -> a -> [a] -> [a]
replace :: forall a. Int -> a -> [a] -> [a]
replace Int
i a
x [a]
xs = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
i [a]
xs [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
xs

swap :: Int -> Int -> [a] -> [a]
swap :: forall a. Int -> Int -> [a] -> [a]
swap Int
i Int
j [a]
xs = Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replace Int
i ([a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
j) (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
replace Int
j ([a]
xs [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) [a]
xs)

epsilon :: (RealFloat a) => a -> a
epsilon :: forall a. RealFloat a => a -> a
epsilon a
x = Integer -> Int -> a
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
1 (Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- a -> Int
forall a. RealFloat a => a -> Int
floatDigits a
x)

defaultTolerance :: (RealFloat a) => Matrix a -> a
defaultTolerance :: forall a. RealFloat a => Matrix a -> a
defaultTolerance Matrix a
a = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
r Int
c) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. RealFloat a => a -> a
epsilon a
prototype
  where
    (Int
r, Int
c) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
a
    prototype :: a
prototype = case Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
a of a
x : [a]
_ -> a
x; [] -> a
0