-- SPDX-License-Identifier: MPL-2.0

{- | Dense real LU factorization with partial row pivoting. This pure reference
implementation accepts 'RealFloat' scalars (including 'Float' and 'Double').
It reuses the existing row-major 'Matrix' API and requires no native backend.

The factors satisfy @P A = L U@. Pivot acceptance uses the largest absolute
entry of the original matrix as a fixed scale. The default relative tolerance
is @n * epsilon@, where epsilon is the gap above one for the scalar type.
This is a pivot-size policy, not a condition-number estimate or an accuracy
guarantee. Use 'luDecomposeWithTolerance' to select a different policy.

All coefficients and right-hand sides must be finite. Arithmetic producing
NaN or infinity returns 'NonFiniteResult'; finite rounding and underflow
follow the scalar type. Existing generic matrix arithmetic is unchanged.
-}
module HMath.Matrix.LU (
    LU,
    luDecompose,
    luDecomposeWithTolerance,
    luLower,
    luUpper,
    luPermutation,
    solveLU,
    solveMatrixLU,
    solve,
    solveMatrix,
    determinantLU,
    determinant,
    inverse,
) where

import Control.Monad (foldM)
import qualified Data.List as List
import HMath.Error (MathError (..))
import HMath.Matrix (
    Matrix,
    fromListChecked,
    identityChecked,
    shape,
    toList,
 )

{- | An immutable, reusable factorization. Its constructor is hidden so solves
can rely on square, finite factors with nonzero pivots and a valid permutation.
-}
data LU a = LU (Matrix a) (Matrix a) [Int] Bool
    deriving (LU a -> LU a -> Bool
(LU a -> LU a -> Bool) -> (LU a -> LU a -> Bool) -> Eq (LU a)
forall a. Eq a => LU a -> LU a -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall a. Eq a => LU a -> LU a -> Bool
== :: LU a -> LU a -> Bool
$c/= :: forall a. Eq a => LU a -> LU a -> Bool
/= :: LU a -> LU a -> Bool
Eq, Int -> LU a -> ShowS
[LU a] -> ShowS
LU a -> String
(Int -> LU a -> ShowS)
-> (LU a -> String) -> ([LU a] -> ShowS) -> Show (LU a)
forall a. Show a => Int -> LU a -> ShowS
forall a. Show a => [LU a] -> ShowS
forall a. Show a => LU a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> LU a -> ShowS
showsPrec :: Int -> LU a -> ShowS
$cshow :: forall a. Show a => LU a -> String
show :: LU a -> String
$cshowList :: forall a. Show a => [LU a] -> ShowS
showList :: [LU a] -> ShowS
Show)

-- | Unit lower-triangular factor.
luLower :: LU a -> Matrix a
luLower :: forall a. LU a -> Matrix a
luLower (LU Matrix a
lower Matrix a
_ [Int]
_ Bool
_) = Matrix a
lower

-- | Upper-triangular factor.
luUpper :: LU a -> Matrix a
luUpper :: forall a. LU a -> Matrix a
luUpper (LU Matrix a
_ Matrix a
upper [Int]
_ Bool
_) = Matrix a
upper

{- | Final zero-based row order: row @i@ of @P A@ is row
@luPermutation factors !! i@ of @A@. This is not a sequence of swap indices.
-}
luPermutation :: LU a -> [Int]
luPermutation :: forall a. LU a -> [Int]
luPermutation (LU Matrix a
_ Matrix a
_ [Int]
permutation Bool
_) = [Int]
permutation

{- | Factor a square matrix using the default relative pivot tolerance.
A @0 x 0@ matrix has a successful empty factorization.
-}
luDecompose :: (RealFloat a) => Matrix a -> Either MathError (LU a)
luDecompose :: forall a. RealFloat a => Matrix a -> Either MathError (LU a)
luDecompose Matrix a
matrix = a -> Matrix a -> Either MathError (LU a)
forall a. RealFloat a => a -> Matrix a -> Either MathError (LU a)
luDecomposeWithTolerance a
tolerance Matrix a
matrix
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
matrix)
    prototype :: a
prototype = case Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
matrix of
        a
value : [a]
_ -> a
value
        [] -> a
0
    tolerance :: a
tolerance = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n a -> a -> a
forall a. Num a => a -> a -> a
* 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
prototype)

{- | Factor with a finite, nonnegative relative tolerance. A zero tolerance
accepts every nonzero pivot. For positive tolerance, a nonzero pivot is
rejected when @abs pivot / maxAbsInput <= tolerance@. Exactly zero pivots
return 'SingularMatrix'; small nonzero pivots return 'NearSingularMatrix'.

Validation order is square shape, tolerance, all coefficient finiteness,
then elimination. Thus a nonfinite entry cannot be hidden by an earlier
singular pivot. Row swaps include the previously computed L multipliers.
-}
luDecomposeWithTolerance ::
    (RealFloat a) => a -> Matrix a -> Either MathError (LU a)
luDecomposeWithTolerance :: forall a. RealFloat a => a -> Matrix a -> Either MathError (LU a)
luDecomposeWithTolerance a
tolerance Matrix a
matrix
    | Int
rows Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
cols = MathError -> Either MathError (LU a)
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
rows, Int
cols))
    | Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
finite a
tolerance) Bool -> Bool -> Bool
|| a
tolerance a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
0 = MathError -> Either MathError (LU a)
forall a b. a -> Either a b
Left MathError
InvalidTolerance
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
validateMatrix Matrix a
matrix
        Int -> [[a]] -> [Int] -> Bool -> Either MathError (LU a)
factor Int
0 (Int -> [a] -> [[a]]
forall a. Int -> [a] -> [[a]]
chunks Int
cols [a]
values) [Int
0 .. Int
rows Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1] Bool
False
  where
    (Int
rows, Int
cols) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
matrix
    values :: [a]
values = Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
matrix
    scale :: a
scale = (a -> a -> a) -> a -> [a] -> a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
List.foldl' (\a
largest a
value -> a -> a -> a
forall a. Ord a => a -> a -> a
max a
largest (a -> a
forall a. Num a => a -> a
abs a
value)) a
0 [a]
values

    factor :: Int -> [[a]] -> [Int] -> Bool -> Either MathError (LU a)
factor Int
k [[a]]
packed [Int]
permutation Bool
oddSwaps
        | Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
rows = do
            lower <- Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
rows Int
rows ([[a]] -> [a]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[a]]
lowerRows)
            upper <- fromListChecked rows rows (concat upperRows)
            Right (LU lower upper permutation oddSwaps)
        | Bool
otherwise = do
            let ([[a]]
finished, [[a]]
candidates) = Int -> [[a]] -> ([[a]], [[a]])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
k [[a]]
packed
                offset :: Int
offset = Int -> [[a]] -> Int
forall a. RealFloat a => Int -> [[a]] -> Int
pivotOffset Int
k [[a]]
candidates
                swapped :: [[a]]
swapped = Int -> Int -> [[a]] -> [[a]]
forall a. Int -> Int -> [a] -> [a]
swapRows Int
0 Int
offset [[a]]
candidates
            case [[a]]
swapped of
                [] -> MathError -> Either MathError (LU a)
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
k)
                [a]
pivotRow : [[a]]
otherRows -> do
                    let pivot :: a
pivot = [a]
pivotRow [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k
                    if a
pivot a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                        then MathError -> Either MathError (LU a)
forall a b. a -> Either a b
Left (Int -> MathError
SingularMatrix Int
k)
                        else
                            if a
tolerance a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 Bool -> Bool -> Bool
&& a -> a
forall a. Num a => a -> a
abs a
pivot a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
scale a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
tolerance
                                then MathError -> Either MathError (LU a)
forall a b. a -> Either a b
Left (Int -> MathError
NearSingularMatrix Int
k)
                                else do
                                    remaining <- ([a] -> Either MathError [a]) -> [[a]] -> Either MathError [[a]]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (Int -> [a] -> [a] -> Either MathError [a]
forall {a}.
RealFloat a =>
Int -> [a] -> [a] -> Either MathError [a]
eliminate Int
k [a]
pivotRow) [[a]]
otherRows
                                    factor
                                        (k + 1)
                                        (finished ++ pivotRow : remaining)
                                        (swapRows k (k + offset) permutation)
                                        (oddSwaps /= (offset /= 0))
      where
        lowerRows :: [[a]]
lowerRows = (Int -> [a] -> [a]) -> [Int] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> [a] -> [a]
forall {a}. Num a => Int -> [a] -> [a]
lowerRow [Int
0 ..] [[a]]
packed
        upperRows :: [[a]]
upperRows = (Int -> [a] -> [a]) -> [Int] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> [a] -> [a]
forall {a}. Num a => Int -> [a] -> [a]
upperRow [Int
0 ..] [[a]]
packed
        lowerRow :: Int -> [a] -> [a]
lowerRow Int
i [a]
row = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
i [a]
row [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> a -> [a]
forall a. Int -> a -> [a]
replicate (Int
rows Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) a
0
        upperRow :: Int -> [a] -> [a]
upperRow Int
i [a]
row = Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
i a
0 [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
i [a]
row

    eliminate :: Int -> [a] -> [a] -> Either MathError [a]
eliminate Int
k [a]
pivotRow [a]
row = do
        multiplier <- a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
checked ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k a -> a -> a
forall a. Fractional a => a -> a -> a
/ ([a]
pivotRow [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
k))
        let update (a
value, a
pivotValue) = do
                productValue <- a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
checked (a
multiplier a -> a -> a
forall a. Num a => a -> a -> a
* a
pivotValue)
                checked (value - productValue)
        trailing <- traverse update (zip (drop (k + 1) row) (drop (k + 1) pivotRow))
        Right (take k row ++ multiplier : trailing)

{- | Solve @A x = b@ using an existing factorization. The RHS must be a finite
list; its length is checked before its values. Each solve preserves the
factorization for reuse.
-}
solveLU :: (RealFloat a) => LU a -> [a] -> Either MathError [a]
solveLU :: forall a. RealFloat a => LU a -> [a] -> Either MathError [a]
solveLU LU a
factors [a]
values
    | Int
actual Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError [a]
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch Int
n Int
actual)
    | Bool
otherwise = do
        [a] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
validateVector [a]
values
        rhs <- Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
n Int
1 [a]
values
        toList <$> solveMatrixLU factors rhs
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape (LU a -> Matrix a
forall a. LU a -> Matrix a
luLower LU a
factors))
    actual :: Int
actual = [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
values

{- | Solve @A X = B@ for any number of right-hand-side columns. Check the RHS
row count before its values. An @n x 0@ RHS and an empty @0 x m@ system retain
their exact shapes, without enumerating empty columns.
-}
solveMatrixLU ::
    (RealFloat a) => LU a -> Matrix a -> Either MathError (Matrix a)
solveMatrixLU :: forall a.
RealFloat a =>
LU a -> Matrix a -> Either MathError (Matrix a)
solveMatrixLU LU a
factors Matrix a
rhs
    | Int
actualRows Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left (Int -> Int -> MathError
RightHandSideMismatch Int
n Int
actualRows)
    | Bool
otherwise = do
        Matrix a -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
validateMatrix Matrix a
rhs
        if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 Bool -> Bool -> Bool
|| Int
columns Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
            then Matrix a -> Either MathError (Matrix a)
forall a b. b -> Either a b
Right Matrix a
rhs
            else do
                intermediate <- [[a]] -> [[a]] -> [[a]] -> Either MathError [[a]]
forall {a}.
RealFloat a =>
[[a]] -> [[a]] -> [[a]] -> Either MathError [[a]]
forward [[a]]
lowerRows [[a]]
permuted []
                solution <- foldM backward [] (reverse (zip upperRows intermediate))
                fromListChecked n columns (concat solution)
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape (LU a -> Matrix a
forall a. LU a -> Matrix a
luLower LU a
factors))
    (Int
actualRows, Int
columns) = Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
rhs
    lowerRows :: [[a]]
lowerRows = Int -> [a] -> [[a]]
forall a. Int -> [a] -> [[a]]
chunks Int
n (Matrix a -> [a]
forall a. Matrix a -> [a]
toList (LU a -> Matrix a
forall a. LU a -> Matrix a
luLower LU a
factors))
    upperRows :: [[a]]
upperRows = Int -> [a] -> [[a]]
forall a. Int -> [a] -> [[a]]
chunks Int
n (Matrix a -> [a]
forall a. Matrix a -> [a]
toList (LU a -> Matrix a
forall a. LU a -> Matrix a
luUpper LU a
factors))
    rhsRows :: [[a]]
rhsRows = Int -> [a] -> [[a]]
forall a. Int -> [a] -> [[a]]
chunks Int
columns (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
rhs)
    permuted :: [[a]]
permuted = (Int -> [a]) -> [Int] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map ([[a]]
rhsRows [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!!) (LU a -> [Int]
forall a. LU a -> [Int]
luPermutation LU a
factors)

    forward :: [[a]] -> [[a]] -> [[a]] -> Either MathError [[a]]
forward [] [] [[a]]
known = [[a]] -> Either MathError [[a]]
forall a b. b -> Either a b
Right [[a]]
known
    forward ([a]
row : [[a]]
remaining) ([a]
b : [[a]]
bs) [[a]]
known = do
        value <- [a] -> [(a, [a])] -> Either MathError [a]
forall a. RealFloat a => [a] -> [(a, [a])] -> Either MathError [a]
subtractRows [a]
b ([a] -> [[a]] -> [(a, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take ([[a]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
known) [a]
row) [[a]]
known)
        forward remaining bs (known ++ [value])
    -- Only internal lists derived from the same validated square order
    -- reach this branch; no user-provided dimensions can reach it.
    forward [[a]]
_ [[a]]
_ [[a]]
_ = MathError -> Either MathError [[a]]
forall a b. a -> Either a b
Left (Int -> Int -> MathError
RightHandSideMismatch Int
n Int
actualRows)

    backward :: [[b]] -> ([b], [b]) -> Either MathError [[b]]
backward [[b]]
known ([b]
row, [b]
b) = do
        let i :: Int
i = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- [[b]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[b]]
known Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
        value <- [b] -> [(b, [b])] -> Either MathError [b]
forall a. RealFloat a => [a] -> [(a, [a])] -> Either MathError [a]
subtractRows [b]
b ([b] -> [[b]] -> [(b, [b])]
forall a b. [a] -> [b] -> [(a, b)]
zip (Int -> [b] -> [b]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [b]
row) [[b]]
known)
        result <- traverse (checked . (/ (row !! i))) value
        Right (result : known)

{- | Factor @A@ with the default tolerance and solve @A x = b@. Matrix
validation and factorization precede RHS validation. For repeated solves or
a custom tolerance, use 'luDecomposeWithTolerance' and 'solveLU'.
-}
solve :: (RealFloat a) => Matrix a -> [a] -> Either MathError [a]
solve :: forall a. RealFloat a => Matrix a -> [a] -> Either MathError [a]
solve Matrix a
matrix [a]
rhs = Matrix a -> Either MathError (LU a)
forall a. RealFloat a => Matrix a -> Either MathError (LU a)
luDecompose Matrix a
matrix Either MathError (LU a)
-> (LU a -> Either MathError [a]) -> Either MathError [a]
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= (LU a -> [a] -> Either MathError [a]
forall a. RealFloat a => LU a -> [a] -> Either MathError [a]
`solveLU` [a]
rhs)

-- | Factor @A@ with the default tolerance and solve @A X = B@.
solveMatrix :: (RealFloat a) => Matrix a -> Matrix a -> Either MathError (Matrix a)
solveMatrix :: forall a.
RealFloat a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
solveMatrix Matrix a
matrix Matrix a
rhs = Matrix a -> Either MathError (LU a)
forall a. RealFloat a => Matrix a -> Either MathError (LU a)
luDecompose Matrix a
matrix Either MathError (LU a)
-> (LU a -> Either MathError (Matrix a))
-> Either MathError (Matrix a)
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= (LU a -> Matrix a -> Either MathError (Matrix a)
forall a.
RealFloat a =>
LU a -> Matrix a -> Either MathError (Matrix a)
`solveMatrixLU` Matrix a
rhs)

{- | Determinant from the permutation sign and diagonal of U. The finite pivots
are multiplied exactly as rationals and rounded once to avoid intermediate
product overflow/underflow. An unrepresentable final result returns
'NonFiniteResult'; final underflow may produce zero and does not imply
singularity. The empty determinant is one.
-}
determinantLU :: (RealFloat a) => LU a -> Either MathError a
determinantLU :: forall a. RealFloat a => LU a -> Either MathError a
determinantLU (LU Matrix a
_ Matrix a
upper [Int]
_ Bool
oddSwaps) = a -> Either MathError a
forall a. RealFloat a => a -> Either MathError a
checked (Rational -> a
forall a. Fractional a => Rational -> a
fromRational Rational
signedProduct)
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
upper)
    diagonal :: [a]
diagonal = ([a] -> Int -> a) -> [[a]] -> [Int] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
(!!) (Int -> [a] -> [[a]]
forall a. Int -> [a] -> [[a]]
chunks Int
n (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
upper)) [Int
0 ..]
    productValue :: Rational
productValue = [Rational] -> Rational
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((a -> Rational) -> [a] -> [Rational]
forall a b. (a -> b) -> [a] -> [b]
map a -> Rational
forall a. Real a => a -> Rational
toRational [a]
diagonal)
    signedProduct :: Rational
signedProduct = if Bool
oddSwaps then Rational -> Rational
forall a. Num a => a -> a
negate Rational
productValue else Rational
productValue

{- | Determinant using the default factorization policy. An exactly zero pivot
returns zero; near-singular and nonfinite failures are propagated. Use an
explicit factorization and 'determinantLU' to choose another tolerance.
-}
determinant :: (RealFloat a) => Matrix a -> Either MathError a
determinant :: forall a. RealFloat a => Matrix a -> Either MathError a
determinant Matrix a
matrix = case Matrix a -> Either MathError (LU a)
forall a. RealFloat a => Matrix a -> Either MathError (LU a)
luDecompose Matrix a
matrix of
    Left (SingularMatrix Int
_) -> a -> Either MathError a
forall a b. b -> Either a b
Right a
0
    Left MathError
err -> MathError -> Either MathError a
forall a b. a -> Either a b
Left MathError
err
    Right LU a
factors -> LU a -> Either MathError a
forall a. RealFloat a => LU a -> Either MathError a
determinantLU LU a
factors

{- | Invert a nonsingular square matrix by solving against the identity. The
inverse of the @0 x 0@ matrix is @0 x 0@. For solving a system directly,
prefer 'solve' or 'solveMatrix' to explicitly forming the inverse.
-}
inverse :: (RealFloat a) => Matrix a -> Either MathError (Matrix a)
inverse :: forall a. RealFloat a => Matrix a -> Either MathError (Matrix a)
inverse Matrix a
matrix = do
    factors <- Matrix a -> Either MathError (LU a)
forall a. RealFloat a => Matrix a -> Either MathError (LU a)
luDecompose Matrix a
matrix
    unit <- identityChecked (fst (shape matrix))
    solveMatrixLU factors unit

subtractRows :: (RealFloat a) => [a] -> [(a, [a])] -> Either MathError [a]
subtractRows :: forall a. RealFloat a => [a] -> [(a, [a])] -> Either MathError [a]
subtractRows = ([a] -> (a, [a]) -> Either MathError [a])
-> [a] -> [(a, [a])] -> Either MathError [a]
forall (t :: * -> *) (m :: * -> *) b a.
(Foldable t, Monad m) =>
(b -> a -> m b) -> b -> t a -> m b
foldM [a] -> (a, [a]) -> Either MathError [a]
forall {b}. RealFloat b => [b] -> (b, [b]) -> Either MathError [b]
subtractRow
  where
    subtractRow :: [b] -> (b, [b]) -> Either MathError [b]
subtractRow [b]
current (b
coefficient, [b]
known) = ((b, b) -> Either MathError b) -> [(b, b)] -> Either MathError [b]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (b, b) -> Either MathError b
update ([b] -> [b] -> [(b, b)]
forall a b. [a] -> [b] -> [(a, b)]
zip [b]
current [b]
known)
      where
        update :: (b, b) -> Either MathError b
update (b
value, b
component) = do
            productValue <- b -> Either MathError b
forall a. RealFloat a => a -> Either MathError a
checked (b
coefficient b -> b -> b
forall a. Num a => a -> a -> a
* b
component)
            checked (value - productValue)

validateMatrix :: (RealFloat a) => Matrix a -> Either MathError ()
validateMatrix :: forall a. RealFloat a => Matrix a -> Either MathError ()
validateMatrix Matrix a
matrix = Int -> [a] -> Either MathError ()
forall {a}. RealFloat a => Int -> [a] -> Either MathError ()
go Int
0 (Matrix a -> [a]
forall a. Matrix a -> [a]
toList Matrix a
matrix)
  where
    columns :: Int
columns = (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix a
matrix)
    go :: Int -> [a] -> Either MathError ()
go Int
_ [] = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    go Int
index (a
value : [a]
rest)
        | a -> Bool
forall a. RealFloat a => a -> Bool
finite a
value = Int -> [a] -> Either MathError ()
go (Int
index Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
rest
        | Bool
otherwise = MathError -> Either MathError ()
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonFiniteMatrixEntry (Int
index Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
columns, Int
index Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
columns))

validateVector :: (RealFloat a) => [a] -> Either MathError ()
validateVector :: forall a. RealFloat a => [a] -> Either MathError ()
validateVector = Int -> [a] -> Either MathError ()
forall {a}. RealFloat a => Int -> [a] -> Either MathError ()
go Int
0
  where
    go :: Int -> [a] -> Either MathError ()
go Int
_ [] = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    go Int
index (a
value : [a]
rest)
        | a -> Bool
forall a. RealFloat a => a -> Bool
finite a
value = Int -> [a] -> Either MathError ()
go (Int
index Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [a]
rest
        | Bool
otherwise = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (Int -> MathError
NonFiniteVectorEntry Int
index)

finite :: (RealFloat a) => a -> Bool
finite :: forall a. RealFloat a => a -> Bool
finite a
value = Bool -> Bool
not (a -> Bool
forall a. RealFloat a => a -> Bool
isNaN a
value Bool -> Bool -> Bool
|| a -> Bool
forall a. RealFloat a => a -> Bool
isInfinite a
value)

checked :: (RealFloat a) => a -> Either MathError a
checked :: forall a. RealFloat a => a -> Either MathError a
checked a
value
    | a -> Bool
forall a. RealFloat a => a -> Bool
finite a
value = a -> Either MathError a
forall a b. b -> Either a b
Right a
value
    | Bool
otherwise = MathError -> Either MathError a
forall a b. a -> Either a b
Left MathError
NonFiniteResult

-- Select the first maximum to make equal-magnitude pivots deterministic.
pivotOffset :: (RealFloat a) => Int -> [[a]] -> Int
pivotOffset :: forall a. RealFloat a => Int -> [[a]] -> Int
pivotOffset Int
_ [] = Int
0
pivotOffset Int
column ([a]
row : [[a]]
rows) = (Int, a) -> Int
forall a b. (a, b) -> a
fst (((Int, a) -> (Int, [a]) -> (Int, a))
-> (Int, a) -> [(Int, [a])] -> (Int, a)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
List.foldl' (Int, a) -> (Int, [a]) -> (Int, a)
forall {b} {a}. (Ord b, Num b) => (a, b) -> (a, [b]) -> (a, b)
choose (Int
0, a -> a
forall a. Num a => a -> a
abs ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
column)) [(Int, [a])]
candidates)
  where
    candidates :: [(Int, [a])]
candidates = [Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 ..] [[a]]
rows
    choose :: (a, b) -> (a, [b]) -> (a, b)
choose best :: (a, b)
best@(a
_, b
largest) (a
index, [b]
candidate)
        | b -> b
forall a. Num a => a -> a
abs ([b]
candidate [b] -> Int -> b
forall a. HasCallStack => [a] -> Int -> a
!! Int
column) b -> b -> Bool
forall a. Ord a => a -> a -> Bool
> b
largest = (a
index, b -> b
forall a. Num a => a -> a
abs ([b]
candidate [b] -> Int -> b
forall a. HasCallStack => [a] -> Int -> a
!! Int
column))
        | Bool
otherwise = (a, b)
best

swapRows :: Int -> Int -> [a] -> [a]
swapRows :: forall a. Int -> Int -> [a] -> [a]
swapRows Int
first Int
second [a]
values
    | Int
first Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
second = [a]
values
    | Bool
otherwise = (Int -> a -> a) -> [Int] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> a -> a
replace [Int
0 ..] [a]
values
  where
    firstValue :: a
firstValue = [a]
values [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
first
    secondValue :: a
secondValue = [a]
values [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
second
    replace :: Int -> a -> a
replace Int
index a
value
        | Int
index Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
first = a
secondValue
        | Int
index Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
second = a
firstValue
        | Bool
otherwise = a
value

chunks :: Int -> [a] -> [[a]]
chunks :: forall a. Int -> [a] -> [[a]]
chunks Int
_ [] = []
chunks Int
width [a]
values = let ([a]
row, [a]
rest) = Int -> [a] -> ([a], [a])
forall a. Int -> [a] -> ([a], [a])
splitAt Int
width [a]
values in [a]
row [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: Int -> [a] -> [[a]]
forall a. Int -> [a] -> [[a]]
chunks Int
width [a]
rest