module Math.Nuha.Algorithms where
import Prelude
import qualified Prelude as P
import Data.Vector.Unboxed (Unbox)
import Math.Nuha.Base
import Math.Nuha.Numeric
import Math.Nuha.Types
solveLin :: (Unbox a, Real a, Floating a)
=> Holor a
-> Holor a -> Either Error (Holor a)
solveLin :: Holor a -> Holor a -> Either Error (Holor a)
solveLin Holor a
_A Holor a
b = case Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isSquare Holor a
_A of
Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
NoSquareMatrixError
Bool
True -> case Holor a -> Holor a -> Either Error (Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Holor a -> Either Error (Holor a)
solveLinLS Holor a
_A Holor a
b of
Left Error
UnderdeterminedSystemError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
UnderdeterminedSystemError
Left Error
DimensionMismatchError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
Right Holor a
x -> Holor a -> Either Error (Holor a)
forall a b. b -> Either a b
Right Holor a
x
solveLinLS :: (Unbox a, Real a, Floating a) => Holor a -> Holor a -> Either Error (Holor a)
solveLinLS :: Holor a -> Holor a -> Either Error (Holor a)
solveLinLS Holor a
_A Holor a
b = case Holor a -> Either Error ([Holor a], Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Either Error ([Holor a], Holor a)
facPreQR Holor a
_A of
Left Error
TooFewRowsError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
UnderdeterminedSystemError
Left Error
NoMatrixError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
Right ([Holor a]
listV, Holor a
_R) -> Either Error (Holor a)
result where
reflection :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
reflection :: Holor a -> Holor a -> Holor a
reflection Holor a
v Holor a
b = Holor a
b Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> a
dot Holor a
v Holor a
b) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| Holor a
v
result :: Either Error (Holor a)
result = case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
b [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
m,Int
1] of
Bool
True -> case Holor a -> Holor a -> Either Error (Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Holor a -> Either Error (Holor a)
solveLinBack Holor a
_R Holor a
_QTb of
Left Error
NoUpperTriError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
UnderdeterminedSystemError
Left Error
DimensionMismatchError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
Left Error
NoMatrixError -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
Right Holor a
x -> Holor a -> Either Error (Holor a)
forall a b. b -> Either a b
Right Holor a
x
Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
where
m :: Int
m = [Int] -> Int
forall a. [a] -> a
head ([Int] -> Int) -> [Int] -> Int
forall a b. (a -> b) -> a -> b
$ Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A
_QTb :: Holor a
_QTb = (Holor a -> Holor a -> Holor a) -> Holor a -> [Holor a] -> Holor a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
reflection Holor a
b [Holor a]
listV
solveLinBack :: (Unbox a, Real a, Floating a) => Holor a -> Holor a -> Either Error (Holor a)
solveLinBack :: Holor a -> Holor a -> Either Error (Holor a)
solveLinBack Holor a
_R Holor a
b = case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_R of
[Int
m,Int
n] -> case Holor a -> Bool
forall a. (Unbox a, Eq a, Num a) => Holor a -> Bool
isUpperTri Holor a
_R of
Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
NoUpperTriError
Bool
True -> case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
b [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
m,Int
1] of
Bool
False -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
Bool
True -> Holor a -> Either Error (Holor a)
forall a b. b -> Either a b
Right (Holor a -> Either Error (Holor a))
-> Holor a -> Either Error (Holor a)
forall a b. (a -> b) -> a -> b
$ Int -> [a] -> Holor a
recursion (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) []
where
recursion :: Int -> [a] -> Holor a
recursion Int
i [a]
xList = case Int
iInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 of
Bool
True -> [a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector [a]
xList'
Bool
False -> Int -> [a] -> Holor a
recursion (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [a]
xList'
where
xList' :: [a]
xList' = a
xi a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
xList
xi :: a
xi = a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/(Holor a
_RHolor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|![Int
i,Int
i]) a -> a -> a
forall a. Num a => a -> a -> a
*
(Holor a
bHolor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|![Int
i,Int
0] a -> a -> a
forall a. Num a => a -> a -> a
-
[a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
P.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
(*) (Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList (Holor a
_RHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||![[Int
i],[Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]])) [a]
xList))
[Int]
_ -> Error -> Either Error (Holor a)
forall a b. a -> Either a b
Left Error
DimensionMismatchError
facQR :: (Unbox a, Real a, Floating a) => Holor a -> Either Error (Holor a, Holor a)
facQR :: Holor a -> Either Error (Holor a, Holor a)
facQR Holor a
_A = case Holor a -> Either Error ([Holor a], Holor a)
forall a.
(Unbox a, Real a, Floating a) =>
Holor a -> Either Error ([Holor a], Holor a)
facPreQR Holor a
_A of
Left Error
TooFewRowsError -> Error -> Either Error (Holor a, Holor a)
forall a b. a -> Either a b
Left Error
TooFewRowsError
Left Error
NoMatrixError -> Error -> Either Error (Holor a, Holor a)
forall a b. a -> Either a b
Left Error
NoMatrixError
Right ([Holor a]
listV, Holor a
_R) -> (Holor a, Holor a) -> Either Error (Holor a, Holor a)
forall a b. b -> Either a b
Right (Holor a
_Q,Holor a
_R) where
_Q :: Holor a
_Q = Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose (Holor a -> Holor a) -> Holor a -> Holor a
forall a b. (a -> b) -> a -> b
$
(Holor a -> Holor a -> Holor a) -> Holor a -> [Holor a] -> Holor a
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\Holor a
v Holor a
_S -> (Int -> Holor a
forall a. (Unbox a, Num a) => Int -> Holor a
eye Int
m Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| (Holor a
v Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|.| Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose Holor a
v)) Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|.| Holor a
_S) (Int -> Holor a
forall a. (Unbox a, Num a) => Int -> Holor a
eye Int
m) [Holor a]
listV
[Int
m,Int
n] = Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A
facPreQR :: (Unbox a, Real a, Floating a) => Holor a -> Either Error ([Holor a], Holor a)
facPreQR :: Holor a -> Either Error ([Holor a], Holor a)
facPreQR Holor a
_A = case Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A of
[Int
m,Int
n] -> Either Error ([Holor a], Holor a)
guards where
guards :: Either Error ([Holor a], Holor a)
guards
| Int
mInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
n = Error -> Either Error ([Holor a], Holor a)
forall a b. a -> Either a b
Left Error
TooFewRowsError
| Int
mInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n = ([Holor a], Holor a) -> Either Error ([Holor a], Holor a)
forall a b. b -> Either a b
Right ([Holor a]
listV, Holor a
_R)
| Bool
otherwise = ([Holor a], Holor a) -> Either Error ([Holor a], Holor a)
forall a b. b -> Either a b
Right ([Holor a]
listV', Holor a
_R')
([Holor a]
listV, [[a]]
listR) = (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
forall a.
(Unbox a, Real a, Floating a) =>
(Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int
m,Int
n) Bool
False [] [] (Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just Holor a
_A)
_R :: Holor a
_R = [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix ([[a]] -> Holor a) -> [[a]] -> Holor a
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [a] -> [a]
reverse [[a]]
listR
([Holor a]
listV', [[a]]
listR') = (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
forall a.
(Unbox a, Real a, Floating a) =>
(Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int
m,Int
n) Bool
True [] [] (Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just Holor a
_A)
_R' :: Holor a
_R' = [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix ([[a]] -> Holor a) -> [[a]] -> Holor a
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [a] -> [a]
reverse (Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
take (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n) ([a] -> [[a]]
forall a. a -> [a]
repeat (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n (a -> [a]
forall a. a -> [a]
repeat a
0))) [[a]] -> [[a]] -> [[a]]
forall a. [a] -> [a] -> [a]
++ [[a]]
listR')
[Int]
_ -> Error -> Either Error ([Holor a], Holor a)
forall a b. a -> Either a b
Left Error
NoMatrixError
where
recursion :: (Unbox a, Real a, Floating a) =>
(Int,Int) -> Bool -> [Holor a] -> [[a]] -> Maybe (Holor a) -> ([Holor a], [[a]])
recursion :: (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int, Int)
shapeIn Bool
additionalStep [Holor a]
listV [[a]]
listR Maybe (Holor a)
maybeA = case Maybe (Holor a)
maybeA of
Maybe (Holor a)
Nothing -> ([Holor a]
listV, [[a]]
listR)
Just Holor a
_A -> (Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
forall a.
(Unbox a, Real a, Floating a) =>
(Int, Int)
-> Bool
-> [Holor a]
-> [[a]]
-> Maybe (Holor a)
-> ([Holor a], [[a]])
recursion (Int, Int)
shapeIn Bool
additionalStep' [Holor a]
listV' [[a]]
listR' Maybe (Holor a)
maybeA' where
(Int
rows,Int
cols) = (Int, Int)
shapeIn
[Int
m,Int
n] = Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
_A
a00 :: a
a00 = Holor a
_AHolor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]
a_0 :: Holor a
a_0 = Holor a
_AHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1],[Int
0]]
sign :: a -> p
sign a
x = if a
xa -> a -> Bool
forall a. Ord a => a -> a -> Bool
<a
0 then -p
1 else p
1
v :: Holor a
v = Holor a -> Holor a
forall a. (Unbox a, Real a, Floating a) => Holor a -> Holor a
normalize2 (Holor a -> Holor a) -> Holor a -> Holor a
forall a b. (a -> b) -> a -> b
$ Holor a
a_0 Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|+| [a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector ((a -> a
forall a p. (Ord a, Num a, Num p) => a -> p
sign a
a00)a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a b. (Unbox a, Real a, Unbox b, Floating b) => Holor a -> b
norm2 Holor a
a_0) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (a -> [a]
forall a. a -> [a]
repeat a
0)))
listV' :: [Holor a]
listV' = ([a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
rows Int -> Int -> Int
forall a. Num a => a -> a -> a
- Holor a -> Int
forall a. Unbox a => Holor a -> Int
numElems Holor a
v) (a -> [a]
forall a. a -> [a]
repeat a
0) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList Holor a
v)) Holor a -> [Holor a] -> [Holor a]
forall a. a -> [a] -> [a]
: [Holor a]
listV
i :: Int
i = [Holor a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Holor a]
listV
([[a]]
listR', Maybe (Holor a)
maybeA', Bool
additionalStep') = case Int
iInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
colsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 of
Bool
True ->
( (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
i (a -> [a]
forall a. a -> [a]
repeat a
0) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ (a
_SA00 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ([[a]] -> [a]
forall a. [a] -> a
head [[a]]
_SA_withoutFirstColumn))) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
listR
, Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just (Holor a -> Maybe (Holor a)) -> Holor a -> Maybe (Holor a)
forall a b. (a -> b) -> a -> b
$ [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix ([[a]] -> [[a]]
forall a. [a] -> [a]
tail [[a]]
_SA_withoutFirstColumn)
, Bool
additionalStep )
Bool
False -> case Bool
additionalStep of
Bool
False ->
( (Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take (Int
colsInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (a -> [a]
forall a. a -> [a]
repeat a
0) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
_SA00]) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
listR
, Maybe (Holor a)
forall a. Maybe a
Nothing
, Bool
False )
Bool
True ->
( [[a]]
listR
, Holor a -> Maybe (Holor a)
forall a. a -> Maybe a
Just (Holor a -> Maybe (Holor a)) -> Holor a -> Maybe (Holor a)
forall a b. (a -> b) -> a -> b
$ Holor a
_SA_firstColumnHolor
, Bool
False )
where
_SA00 :: a
_SA00 = - (a -> a
forall a p. (Ord a, Num a, Num p) => a -> p
sign a
a00) a -> a -> a
forall a. Num a => a -> a -> a
* (Holor a -> a
forall a b. (Unbox a, Real a, Unbox b, Floating b) => Holor a -> b
norm2 Holor a
a_0)
_SA_withoutFirstColumn :: [[a]]
_SA_withoutFirstColumn = Holor a -> [[a]]
forall a. Unbox a => Holor a -> [[a]]
toList2 (Holor a -> [[a]]) -> Holor a -> [[a]]
forall a b. (a -> b) -> a -> b
$ Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose (Holor a -> Holor a) -> Holor a -> Holor a
forall a b. (a -> b) -> a -> b
$ [[a]] -> Holor a
forall a. Unbox a => [[a]] -> Holor a
matrix [
let
a :: Holor a
a = Holor a
_AHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1],[Int
j]]
in
Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList (Holor a -> [a]) -> Holor a -> [a]
forall a b. (a -> b) -> a -> b
$ Holor a
a Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> a
dot Holor a
v Holor a
a) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| Holor a
v
| Int
j <- [Int
1..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]]
_SA_firstColumnHolor :: Holor a
_SA_firstColumnHolor = Holor a
a Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|-| a
2a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> a
dot Holor a
v Holor a
a) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| Holor a
v where
a :: Holor a
a = Holor a
_AHolor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1],[Int
0]]