-- {-# LANGUAGE BangPatterns #-} -- for Debug.Trace

-- |
-- Copyright   : (c) Johannes Kropp
-- License     : BSD 3-Clause
-- Maintainer  : Johannes Kropp <jodak932@gmail.com>

module Math.Nuha.Numeric where

import qualified Prelude as P
import Prelude hiding (map, sum, replicate)
import Data.Vector.Unboxed (Unbox)
import qualified Data.Vector.Unboxed as V
-- import qualified Debug.Trace as D

import Math.Nuha.Types
import Math.Nuha.Base
import Math.Nuha.Internal


{- TODO:
Is there room for optimize things?. See:
    Paper: "Exploiting Vector Instructions with Generalized Stream Fusion"
    https://gitlab.haskell.org/ghc/ghc/-/wikis/simd
    https://gitlab.haskell.org/ghc/ghc/-/wikis/simd/design
    https://github.com/haskell/vector/issues/251
    https://github.com/Magalame/lascive
-}

----------------
-- ** Constructions

-- | diagonalizes a list of values as a square matrix
diag :: (Unbox a, Num a) => [a] -> Holor a
{-# INLINE diag #-}
diag :: [a] -> Holor a
diag [a]
lst = [[Int]] -> [a] -> Holor a -> Holor a
forall a. Unbox a => [[Int]] -> [a] -> Holor a -> Holor a
setElems [[Int]]
mIdcs [a]
lst Holor a
hlr0 where
    len :: Int
len = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
lst
    shape :: [Int]
shape = [Int
len, Int
len]
    hlr0 :: Holor a
hlr0 = [Int] -> Holor a
forall a. (Unbox a, Num a) => [Int] -> Holor a
zeros [Int
len, Int
len]
    mIdcs :: [[Int]]
mIdcs = [[Int
i,Int
i] | Int
i <- [Int
0..Int
lenInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]]

-- | holor with all entries zero
zeros :: (Unbox a, Num a) => [Int] -> Holor a
{-# INLINE zeros #-}
zeros :: [Int] -> Holor a
zeros [Int]
shape = [Int] -> a -> Holor a
forall a. Unbox a => [Int] -> a -> Holor a
replicate [Int]
shape a
0

-- | holor with all entries one
ones :: (Unbox a, Num a) => [Int] -> Holor a
{-# INLINE ones #-}
ones :: [Int] -> Holor a
ones [Int]
shape = [Int] -> a -> Holor a
forall a. Unbox a => [Int] -> a -> Holor a
replicate [Int]
shape a
1

-- | identity matrix
eye :: (Unbox a, Num a) => Int -> Holor a
{-# INLINE eye #-}
eye :: Int -> Holor a
eye Int
dim = [a] -> Holor a
forall a. (Unbox a, Num a) => [a] -> Holor a
diag ([a] -> Holor a) -> [a] -> Holor a
forall a b. (a -> b) -> a -> b
$ Int -> a -> [a]
forall a. Int -> a -> [a]
P.replicate Int
dim a
1

{- | creates a column vector with points in linear sequence

>>> linSpace 1 3 5
  1.0
  1.5
  2.0
  2.5
  3.0
-}
linSpace :: (Unbox a, Fractional a, Enum a) => a -> a -> Int -> Holor a
{-# INLINE linSpace #-}
linSpace :: a -> a -> Int -> Holor a
linSpace a
start a
stop Int
num = [a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector [a]
values where
    values :: [a]
values = [a
start a -> a -> a
forall a. Num a => a -> a -> a
+ a
ia -> a -> a
forall a. Num a => a -> a -> a
*a
step | a
i <- [a
0.0 .. Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral(Int
numInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)]]
    step :: a
step = (a
stop a -> a -> a
forall a. Num a => a -> a -> a
- a
start)a -> a -> a
forall a. Fractional a => a -> a -> a
/(Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
numInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))


---------------
-- ** Operators

infixl 6 |+
infixl 6 +|
infixl 6 |-
infixl 6 -|
infixr 7 |*
infixl 7 *|
infixr 7 |/
infixr 7 /|
infixl 6 |+|
infixl 6 |-|
infixl 7 |*|
infixl 7 |.|
infixr 8 |**
infixr 8 |^

-- | Add a scalar from the right to each element of a holor
(|+) :: (Unbox a, Num a) => Holor a -> a -> Holor a
{-# INLINE (|+) #-}
|+ :: Holor a -> a -> Holor a
(|+) Holor a
hlr a
scalar = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Num a => a -> a -> a
+a
scalar) Holor a
hlr

-- | Add a scalar from the left to each element of a holor
(+|) :: (Unbox a, Num a) => a -> Holor a -> Holor a
{-# INLINE (+|) #-}
+| :: a -> Holor a -> Holor a
(+|) a
scalar Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Num a => a -> a -> a
+a
scalar) Holor a
hlr

-- | Subtract a scalar from the right to each element of a holor
(|-) :: (Unbox a, Num a) => Holor a -> a -> Holor a
{-# INLINE (|-) #-}
|- :: Holor a -> a -> Holor a
(|-) Holor a
hlr a
scalar = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Num a => a -> a -> a
+(-a
scalar)) Holor a
hlr

-- | A single scalar is subtracted by each element of a holor
(-|) :: (Unbox a, Num a) => a -> Holor a -> Holor a
{-# INLINE (-|) #-}
-| :: a -> Holor a -> Holor a
(-|) a
scalar Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a
scalara -> a -> a
forall a. Num a => a -> a -> a
-) Holor a
hlr

-- | Holor multiplied by a scalar from right
(|*) :: (Unbox a, Num a) => Holor a -> a -> Holor a
{-# INLINE (|*) #-}
|* :: Holor a -> a -> Holor a
(|*) Holor a
hlr a
factor = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Num a => a -> a -> a
*a
factor) Holor a
hlr

-- | Holor multiplied by a scalar from left
(*|) :: (Unbox a, Num a) => a -> Holor a -> Holor a
{-# INLINE (*|) #-}
*| :: a -> Holor a -> Holor a
(*|) a
factor Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Num a => a -> a -> a
*a
factor) Holor a
hlr

-- | Holor divided by a scalar
(|/) :: (Unbox a, Fractional a) => Holor a -> a -> Holor a
{-# INLINE (|/) #-}
|/ :: Holor a -> a -> Holor a
(|/) Holor a
hlr a
divisor = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/a
divisor) Holor a
hlr

-- | A single scalar is divided by each element of a holor
(/|) :: (Unbox a, Fractional a) => a -> Holor a -> Holor a
{-# INLINE (/|) #-}
/| :: a -> Holor a -> Holor a
(/|) a
dividend Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a
dividenda -> a -> a
forall a. Fractional a => a -> a -> a
/) Holor a
hlr

-- | elementwise holor addition
(|+|) :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
{-# INLINE (|+|) #-}
|+| :: Holor a -> Holor a -> Holor a
(|+|) Holor a
h1 Holor a
h2
    | (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h1) [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h2) =
        [Int] -> [Int] -> Vector a -> Holor a
forall a. [Int] -> [Int] -> Vector a -> Holor a
Holor (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h1) (Holor a -> [Int]
forall a. Holor a -> [Int]
hStrides Holor a
h1) ((a -> a -> a) -> Vector a -> Vector a -> Vector a
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h1) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h2))
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|+|) : Shape missmatch at elementwise holor addition"

-- | elementwise holor subtraction
(|-|) :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
{-# INLINE (|-|) #-}
|-| :: Holor a -> Holor a -> Holor a
(|-|) Holor a
h1 Holor a
h2
    | (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h1) [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h2) =
        [Int] -> [Int] -> Vector a -> Holor a
forall a. [Int] -> [Int] -> Vector a -> Holor a
Holor (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h1) (Holor a -> [Int]
forall a. Holor a -> [Int]
hStrides Holor a
h1) ((a -> a -> a) -> Vector a -> Vector a -> Vector a
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (-) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h1) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h2))
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|-|) : Shape missmatch at elementwise holor subtraction"

-- | elementwise holor multiplication
(|*|) :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
{-# INLINE (|*|) #-}
|*| :: Holor a -> Holor a -> Holor a
(|*|) Holor a
h1 Holor a
h2
    | (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h1) [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h2) =
        [Int] -> [Int] -> Vector a -> Holor a
forall a. [Int] -> [Int] -> Vector a -> Holor a
Holor (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
h1) (Holor a -> [Int]
forall a. Holor a -> [Int]
hStrides Holor a
h1) ((a -> a -> a) -> Vector a -> Vector a -> Vector a
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h1) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h2))
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|*|) : Shape missmatch at elementwise holor multiplication"


-- | Holor multiplication, last dimension of left holor has to match first dimension of right holor
(|.|) :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
{-# INLINE (|.|) #-}
|.| :: Holor a -> Holor a -> Holor a
(|.|) Holor a
h1 Holor a
h2
    | Holor a -> Int
forall a. Unbox a => Holor a -> Int
dim Holor a
h1 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 Bool -> Bool -> Bool
|| Holor a -> Int
forall a. Unbox a => Holor a -> Int
dim Holor a
h2 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|.|) : holor dimensions have to be at least 2"
    | [Int] -> Int
forall a. [a] -> a
last (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h1) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
matchdim =
        [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|.|) : dimension mismatch: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h1) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" , " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h2)
    | Bool
otherwise = [Int] -> [Int] -> Vector a -> Holor a
forall a. [Int] -> [Int] -> Vector a -> Holor a
Holor [Int]
shapeOut [Int]
stridesOut Vector a
valuesOut
    where
        shape1 :: [Int]
shape1 = Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
take ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h1) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h1)
        -- !shape1_ = D.trace ("shape1: " ++ show shape1) ()
        shape2 :: [Int]
shape2 = Int -> [Int] -> [Int]
forall a. Int -> [a] -> [a]
drop Int
1 (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h2)
        -- !shape2_ = D.trace ("shape2: " ++ show shape2) ()
        shapeOut :: [Int]
shapeOut = [Int]
shape1 [Int] -> [Int] -> [Int]
forall a. [a] -> [a] -> [a]
++ [Int]
shape2
        stridesOut :: [Int]
stridesOut = [Int] -> [Int]
fromShapeToStrides [Int]
shapeOut
        matchdim :: Int
matchdim = [Int] -> Int
forall a. [a] -> a
head (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
h2) -- also step
        -- !matchdim_ = D.trace ("matchdim: " ++ show matchdim) ()
        h2T :: Holor a
h2T = Holor a -> Holor a
forall a. Unbox a => Holor a -> Holor a
transpose Holor a
h2
        values1 :: Vector a
values1 = Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h1
        values2 :: Vector a
values2 = Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
h2T
        valuesOut :: Vector a
valuesOut = [a] -> Vector a
forall a. Unbox a => [a] -> Vector a
V.fromList [
            let
                slice1 :: Vector a
slice1 = Int -> Int -> Vector a -> Vector a
forall a. Unbox a => Int -> Int -> Vector a -> Vector a
V.unsafeSlice (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
matchdim) Int
matchdim Vector a
values1
                slice2 :: Vector a
slice2 = Int -> Int -> Vector a -> Vector a
forall a. Unbox a => Int -> Int -> Vector a -> Vector a
V.unsafeSlice (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
matchdim) Int
matchdim Vector a
values2
            in
                Vector a -> a
forall a. (Unbox a, Num a) => Vector a -> a
V.sum ((a -> a -> a) -> Vector a -> Vector a -> Vector a
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) Vector a
slice1 Vector a
slice2)
            | Int
i <- [Int
0 .. ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int]
shape1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)], Int
j <- [Int
0 .. ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int]
shape2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)]]


-- | Holor to a power elementwise
(|**) :: (Unbox a, Floating a) => Holor a -> a -> Holor a
{-# INLINE (|**) #-}
|** :: Holor a -> a -> Holor a
(|**) Holor a
hlr a
p = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Floating a => a -> a -> a
**a
p) Holor a
hlr

-- TODO: may be optimized by saving power of two products
-- | Holor multiplied with itself elementwise n-times
(|^) :: (Unbox a, Num a) => Holor a -> Int -> Holor a
{-# INLINE (|^) #-}
|^ :: Holor a -> Int -> Holor a
(|^) Holor a
hlr Int
num
    | Int
numInt -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<Int
0 = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|^) : number of multiplications must be non negative"
    | Bool -> Bool
not (Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isSquare Holor a
hlr) = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"(|^) : holor is not of square shape :" [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show (Holor a -> [Int]
forall a. Holor a -> [Int]
hShape Holor a
hlr)
    | Int
num Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = Int -> Holor a
forall a. (Unbox a, Num a) => Int -> Holor a
eye (Int -> Holor a) -> Int -> Holor a
forall a b. (a -> b) -> a -> b
$ (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
hlr) [Int] -> Int -> Int
forall a. [a] -> Int -> a
P.!! Int
0
    | Bool
otherwise = Holor a -> Int -> Holor a
iterMult Holor a
hlr Int
0
    where
        iterMult :: Holor a -> Int -> Holor a
iterMult Holor a
hlr_ Int
n
            | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = Holor a -> Int -> Holor a
iterMult Holor a
hlr (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
            | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
num = Holor a -> Int -> Holor a
iterMult (Holor a
hlr Holor a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a -> Holor a
|.| Holor a
hlr_) (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
            | Bool
otherwise = Holor a
hlr_

----------------
-- ** Operations

-- | additive inverse of a holor
negate :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE negate #-}
negate :: Holor a -> Holor a
negate = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map a -> a
forall a. Num a => a -> a
P.negate

-- | absolute values ​​by holor elements
abs :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE abs #-}
abs :: Holor a -> Holor a
abs= (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map a -> a
forall a. Num a => a -> a
P.abs

-- | signum by holor elements
signum :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE signum #-}
signum :: Holor a -> Holor a
signum = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map a -> a
forall a. Num a => a -> a
P.signum

-- | maximum of all holor elements
maximum :: (Unbox a, Ord a) => Holor a -> a
{-# INLINE maximum #-}
maximum :: Holor a -> a
maximum Holor a
hlr = Vector a -> a
forall a. (Unbox a, Ord a) => Vector a -> a
V.maximum (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr

-- | minimum of all holor elements
minimum :: (Unbox a, Ord a) => Holor a -> a
{-# INLINE minimum #-}
minimum :: Holor a -> a
minimum Holor a
hlr = Vector a -> a
forall a. (Unbox a, Ord a) => Vector a -> a
V.minimum (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr

-- | sum all element of a holor
sum :: (Unbox a, Num a) => Holor a -> a
{-# INLINE sum #-}
sum :: Holor a -> a
sum Holor a
hlr = Vector a -> a
forall a. (Unbox a, Num a) => Vector a -> a
V.sum (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr

-- | trace of a matrix
trace :: (Unbox a, Num a) => Holor a -> a
{-# INLINE trace #-}
trace :: Holor a -> a
trace Holor a
hlr
    | Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isSquare Holor a
hlr = Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
sum (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
diagonal Holor a
hlr
    | Bool
otherwise = [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char] -> a) -> [Char] -> a
forall a b. (a -> b) -> a -> b
$
        [Char]
"trace : holor with shape " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show (Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
hlr) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" is not a square matrix"

-- | p-norm of a holor
norm :: (Unbox a, Real a, Unbox b, Floating b) => a -> Holor a -> b
{-# INLINE norm #-}
norm :: a -> Holor a -> b
norm a
p Holor a
hlr = (Vector b -> b
forall a. (Unbox a, Num a) => Vector a -> a
V.sum (Vector b -> b) -> Vector b -> b
forall a b. (a -> b) -> a -> b
$ (a -> b) -> Vector a -> Vector b
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\a
v -> (a -> b
forall a b. (Real a, Fractional b) => a -> b
realToFrac a
v)b -> b -> b
forall a. Floating a => a -> a -> a
**b
p') (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr))b -> b -> b
forall a. Floating a => a -> a -> a
**(b
1b -> b -> b
forall a. Fractional a => a -> a -> a
/b
p') where
    p' :: b
p' = a -> b
forall a b. (Real a, Fractional b) => a -> b
realToFrac a
p

-- | manhattan norm of a holor (p = 1)
norm1 :: (Unbox a, Num a) => Holor a -> a
{-# INLINE norm1 #-}
norm1 :: Holor a -> a
norm1 Holor a
hlr = Vector a -> a
forall a. (Unbox a, Num a) => Vector a -> a
V.sum (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a) -> Vector a -> Vector a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\a
v -> a -> a
forall a. Num a => a -> a
P.abs a
v) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr)

-- | frobenius norm of a holor. If the holor has vector form it's equal to the euclidian norm
norm2 :: (Unbox a, Real a, Unbox b, Floating b) => Holor a -> b
{-# INLINE norm2 #-}
norm2 :: Holor a -> b
norm2 Holor a
hlr = b -> b
forall a. Floating a => a -> a
sqrt (b -> b) -> b -> b
forall a b. (a -> b) -> a -> b
$ Vector b -> b
forall a. (Unbox a, Num a) => Vector a -> a
V.sum (Vector b -> b) -> Vector b -> b
forall a b. (a -> b) -> a -> b
$ (a -> b) -> Vector a -> Vector b
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\a
v -> (a -> b
forall a b. (Real a, Fractional b) => a -> b
realToFrac a
v)b -> b -> b
forall a. Floating a => a -> a -> a
**b
2) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr)

-- | maximum norm of a holor (p = inf)
normM :: (Unbox a, Num a, Ord a) => Holor a -> a
{-# INLINE normM #-}
normM :: Holor a -> a
normM Holor a
hlr = Vector a -> a
forall a. (Unbox a, Ord a) => Vector a -> a
V.maximum (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a) -> Vector a -> Vector a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\a
v -> a -> a
forall a. Num a => a -> a
P.abs a
v) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
hlr)

-- | normalize a holor to unit-length in p-norm
normalize :: (Unbox a, Real a, Floating a) => a -> Holor a -> Holor a
{-# INLINE normalize #-}
normalize :: a -> Holor a -> Holor a
normalize a
p Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/(a -> Holor a -> a
forall a b.
(Unbox a, Real a, Unbox b, Floating b) =>
a -> Holor a -> b
norm a
p Holor a
hlr)) Holor a
hlr

-- | normalize a holor to unit-length in 1-norm
normalize1 :: (Unbox a, Fractional a) => Holor a -> Holor a
{-# INLINE normalize1 #-}
normalize1 :: Holor a -> Holor a
normalize1 Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
norm1 Holor a
hlr)) Holor a
hlr

-- | normalize a holor to unit-length in 2-norm
normalize2 :: (Unbox a, Real a, Floating a) => Holor a -> Holor a
{-# INLINE normalize2 #-}
normalize2 :: Holor a -> Holor a
normalize2 Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/(Holor a -> a
forall a b. (Unbox a, Real a, Unbox b, Floating b) => Holor a -> b
norm2 Holor a
hlr)) Holor a
hlr
-- | normalize a holor to unit-length in maximum-norm

normalizeM :: (Unbox a, Fractional a, Ord a) => Holor a -> Holor a
{-# INLINE normalizeM #-}
normalizeM :: Holor a -> Holor a
normalizeM Holor a
hlr = (a -> a) -> Holor a -> Holor a
forall a b. (Unbox a, Unbox b) => (a -> b) -> Holor a -> Holor b
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/(Holor a -> a
forall a. (Unbox a, Num a, Ord a) => Holor a -> a
normM Holor a
hlr)) Holor a
hlr

-- | dot product for holors with shape [n,1] (column-vector) or [1,n] (row-vector)
dot :: (Unbox a, Num a) => Holor a -> Holor a -> a
{-# INLINE dot #-}
dot :: Holor a -> Holor a -> a
dot Holor a
v1 Holor a
v2
    | Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isVector Holor a
v1 Bool -> Bool -> Bool
&& Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isVector Holor a
v2 =
        [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char] -> a) -> [Char] -> a
forall a b. (a -> b) -> a -> b
$ [Char]
"dot : arguments have no vector form"
    | [Int]
shape1 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
/= [Int]
shape2 =
        [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char] -> a) -> [Char] -> a
forall a b. (a -> b) -> a -> b
$ [Char]
"dot : shapes " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show ([Int]
shape1) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" and " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show ([Int]
shape2) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" incompatible"
    | Bool
otherwise =
        Vector a -> a
forall a. (Unbox a, Num a) => Vector a -> a
V.sum (Vector a -> a) -> Vector a -> a
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> Vector a -> Vector a -> Vector a
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
v1) (Holor a -> Vector a
forall a. Holor a -> Vector a
hValues Holor a
v2)
    where
        shape1 :: [Int]
shape1 = Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
v1
        shape2 :: [Int]
shape2 = Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
v2

-- | cross product for holors with shape [3,1] (column-vector) or [1,3] (row-vector)
cross :: (Unbox a, Num a) => Holor a -> Holor a -> Holor a
{-# INLINE cross #-}
cross :: Holor a -> Holor a -> Holor a
cross Holor a
v1 Holor a
v2
    | Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isVector Holor a
v1 Bool -> Bool -> Bool
&& Holor a -> Bool
forall a. Unbox a => Holor a -> Bool
isVector Holor a
v2 =
        [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"cross : both arguments must be a column or row vector"
    | Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ [Int]
shape1 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int]
shape2 Bool -> Bool -> Bool
&& ([Int]
shape1 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
3,Int
1] Bool -> Bool -> Bool
|| [Int]
shape1 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
1,Int
3]) =
        [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char] -> Holor a) -> [Char] -> Holor a
forall a b. (a -> b) -> a -> b
$ [Char]
"cross : shapes " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show ([Int]
shape1) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" and " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Int] -> [Char]
forall a. Show a => a -> [Char]
show ([Int]
shape2) [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
" incompatible"
    | Bool
otherwise = [a] -> Holor a
forall a. Unbox a => [a] -> Holor a
vector [a
v11a -> a -> a
forall a. Num a => a -> a -> a
*a
v22a -> a -> a
forall a. Num a => a -> a -> a
-a
v12a -> a -> a
forall a. Num a => a -> a -> a
*a
v21, a
v12a -> a -> a
forall a. Num a => a -> a -> a
*a
v20a -> a -> a
forall a. Num a => a -> a -> a
-a
v10a -> a -> a
forall a. Num a => a -> a -> a
*a
v22, a
v10a -> a -> a
forall a. Num a => a -> a -> a
*a
v21a -> a -> a
forall a. Num a => a -> a -> a
-a
v11a -> a -> a
forall a. Num a => a -> a -> a
*a
v20]
    where
        shape1 :: [Int]
shape1 = Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
v1
        shape2 :: [Int]
shape2 = Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
v2
        [a
v10, a
v11, a
v12] = Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList Holor a
v1
        [a
v20, a
v21, a
v22] = Holor a -> [a]
forall a. Unbox a => Holor a -> [a]
toList Holor a
v2

{- TODO:
-- Link: https://en.wikipedia.org/wiki/Adjugate_matrix
-- | adjugate Matrix
adjugate :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE adjugate #-}
adjugate hlr = undefined
-}

-- | adjugate of 2x2 matrix
adjugate22 :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE adjugate22 #-}
adjugate22 :: Holor a -> Holor a
adjugate22 Holor a
mat22
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat22 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
2,Int
2] = [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
2,Int
2] [a]
values
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char]
"adjugate22 only possible for a holor with shape [2,2]")
    where
        m00 :: a
m00 =   Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]
        m01 :: a
m01 = - Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
0]
        m10 :: a
m10 = - Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
1]
        m11 :: a
m11 =   Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
1]
        values :: [a]
values = [a
m11, a
m10, a
m01, a
m00]

-- | adjugate of 3x3 matrix
adjugate33 :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE adjugate33 #-}
adjugate33 :: Holor a -> Holor a
adjugate33 Holor a
mat33
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat33 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
3,Int
3] = [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a]
values
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char]
"adjugate33 only possible for a holor with shape [3,3]")
    where
        h00 :: a
h00 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]; h01 :: a
h01 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
1]; h02 :: a
h02 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
2]
        h10 :: a
h10 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
0]; h11 :: a
h11 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
1]; h12 :: a
h12 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
2]
        h20 :: a
h20 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
0]; h21 :: a
h21 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
1]; h22 :: a
h22 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
2]
        m00 :: a
m00 = a
h11a -> a -> a
forall a. Num a => a -> a -> a
*a
h22 a -> a -> a
forall a. Num a => a -> a -> a
- a
h12a -> a -> a
forall a. Num a => a -> a -> a
*a
h21
        m01 :: a
m01 = a
h12a -> a -> a
forall a. Num a => a -> a -> a
*a
h20 a -> a -> a
forall a. Num a => a -> a -> a
- a
h10a -> a -> a
forall a. Num a => a -> a -> a
*a
h22
        m02 :: a
m02 = a
h10a -> a -> a
forall a. Num a => a -> a -> a
*a
h21 a -> a -> a
forall a. Num a => a -> a -> a
- a
h11a -> a -> a
forall a. Num a => a -> a -> a
*a
h20
        m10 :: a
m10 = a
h02a -> a -> a
forall a. Num a => a -> a -> a
*a
h21 a -> a -> a
forall a. Num a => a -> a -> a
- a
h01a -> a -> a
forall a. Num a => a -> a -> a
*a
h22
        m11 :: a
m11 = a
h00a -> a -> a
forall a. Num a => a -> a -> a
*a
h22 a -> a -> a
forall a. Num a => a -> a -> a
- a
h02a -> a -> a
forall a. Num a => a -> a -> a
*a
h20
        m12 :: a
m12 = a
h01a -> a -> a
forall a. Num a => a -> a -> a
*a
h20 a -> a -> a
forall a. Num a => a -> a -> a
- a
h00a -> a -> a
forall a. Num a => a -> a -> a
*a
h21
        m20 :: a
m20 = a
h01a -> a -> a
forall a. Num a => a -> a -> a
*a
h12 a -> a -> a
forall a. Num a => a -> a -> a
- a
h02a -> a -> a
forall a. Num a => a -> a -> a
*a
h11
        m21 :: a
m21 = a
h02a -> a -> a
forall a. Num a => a -> a -> a
*a
h10 a -> a -> a
forall a. Num a => a -> a -> a
- a
h00a -> a -> a
forall a. Num a => a -> a -> a
*a
h12
        m22 :: a
m22 = a
h00a -> a -> a
forall a. Num a => a -> a -> a
*a
h11 a -> a -> a
forall a. Num a => a -> a -> a
- a
h01a -> a -> a
forall a. Num a => a -> a -> a
*a
h10
        values :: [a]
values = [a
m00, a
m10, a
m20, a
m01, a
m11, a
m21, a
m02, a
m12, a
m22]

-- Link: https://semath.info/src/inverse-cofactor-ex4.html
-- (Adj44)_ij = (-1)^(i+j) * det33(M_ji)
-- | adjugate of 4x4 matrix
adjugate44 :: (Unbox a, Num a) => Holor a -> Holor a
{-# INLINE adjugate44 #-}
adjugate44 :: Holor a -> Holor a
adjugate44 Holor a
mat44
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat44 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
4,Int
4] = [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
4,Int
4] [a]
values
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char]
"adjugate44 only possible for a holor with shape [4,4]")
    where
        h00 :: a
h00 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]; h01 :: a
h01 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
1]; h02 :: a
h02 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
2]; h03 :: a
h03 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
3]
        h10 :: a
h10 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
0]; h11 :: a
h11 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
1]; h12 :: a
h12 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
2]; h13 :: a
h13 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
3]
        h20 :: a
h20 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
0]; h21 :: a
h21 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
1]; h22 :: a
h22 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
2]; h23 :: a
h23 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
2,Int
3]
        h30 :: a
h30 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
3,Int
0]; h31 :: a
h31 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
3,Int
1]; h32 :: a
h32 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
3,Int
2]; h33 :: a
h33 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
3,Int
3]
        m00 :: a
m00 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h11,a
h12,a
h13,a
h21,a
h22,a
h23,a
h31,a
h32,a
h33])
        m01 :: a
m01 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h10,a
h12,a
h13,a
h20,a
h22,a
h23,a
h30,a
h32,a
h33])
        m02 :: a
m02 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h10,a
h11,a
h13,a
h20,a
h21,a
h23,a
h30,a
h31,a
h33])
        m03 :: a
m03 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h10,a
h11,a
h12,a
h20,a
h21,a
h22,a
h30,a
h31,a
h32])
        m10 :: a
m10 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h01,a
h02,a
h03,a
h21,a
h22,a
h23,a
h31,a
h32,a
h33])
        m11 :: a
m11 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h02,a
h03,a
h20,a
h22,a
h23,a
h30,a
h32,a
h33])
        m12 :: a
m12 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h01,a
h03,a
h20,a
h21,a
h23,a
h30,a
h31,a
h33])
        m13 :: a
m13 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h01,a
h02,a
h20,a
h21,a
h22,a
h30,a
h31,a
h32])
        m20 :: a
m20 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h01,a
h02,a
h03,a
h11,a
h12,a
h13,a
h31,a
h32,a
h33])
        m21 :: a
m21 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h02,a
h03,a
h10,a
h12,a
h13,a
h30,a
h32,a
h33])
        m22 :: a
m22 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h01,a
h03,a
h10,a
h11,a
h13,a
h30,a
h31,a
h33])
        m23 :: a
m23 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h01,a
h02,a
h10,a
h11,a
h12,a
h30,a
h31,a
h32])
        m30 :: a
m30 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h01,a
h02,a
h03,a
h11,a
h12,a
h13,a
h21,a
h22,a
h23])
        m31 :: a
m31 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h02,a
h03,a
h10,a
h12,a
h13,a
h20,a
h22,a
h23])
        m32 :: a
m32 = - (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h01,a
h03,a
h10,a
h11,a
h13,a
h20,a
h21,a
h23])
        m33 :: a
m33 =   (Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 (Holor a -> a) -> Holor a -> a
forall a b. (a -> b) -> a -> b
$ [Int] -> [a] -> Holor a
forall a. Unbox a => [Int] -> [a] -> Holor a
holor [Int
3,Int
3] [a
h00,a
h01,a
h02,a
h10,a
h11,a
h12,a
h20,a
h21,a
h22])
        values :: [a]
values = [a
m00,a
m10,a
m20,a
m30,a
m01,a
m11,a
m21,a
m31,a
m02,a
m12,a
m22,a
m32,a
m03,a
m13,a
m23,a
m33]

{- TODO:
det :: (Unbox a, Num a) => Holor a -> a
{-# INLINE det #-}
det hlr
    | (length (shape hlr) == 2) && ((shape hlr) P.!! 0 == (shape hlr) P.!! 1 = result
    | otherwise = error $ "det of nonsquare holor with shape " ++ show (shape hlr)
    where
        result = undefined
-}


-- | determinant of 2x2 matrix
det22 :: (Unbox a, Num a) => Holor a -> a
{-# INLINE det22 #-}
det22 :: Holor a -> a
det22 Holor a
mat22
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat22 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
2,Int
2] = a
h00a -> a -> a
forall a. Num a => a -> a -> a
*a
h11 a -> a -> a
forall a. Num a => a -> a -> a
- a
h01a -> a -> a
forall a. Num a => a -> a -> a
*a
h10
    | Bool
otherwise = [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char]
"det22 only possible for a holor with shape [2,2]")
    where
        h00 :: a
h00 = Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]
        h01 :: a
h01 = Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
1]
        h10 :: a
h10 = Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
0]
        h11 :: a
h11 = Holor a
mat22Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
1,Int
1]

-- Laplace-Expansion by first row
-- | determinant of 3x3 matrix
det33 :: (Unbox a, Num a) => Holor a -> a
{-# INLINE det33 #-}
det33 :: Holor a -> a
det33 Holor a
mat33
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat33 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
3,Int
3] = a
h00a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det22 Holor a
m00) a -> a -> a
forall a. Num a => a -> a -> a
- a
h01a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det22 Holor a
m01) a -> a -> a
forall a. Num a => a -> a -> a
+ a
h02a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det22 Holor a
m02)
    | Bool
otherwise = [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char]
"det33 only possible for a holor with shape [3,3]")
    where
        h00 :: a
h00 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]
        h01 :: a
h01 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
1]
        h02 :: a
h02 = Holor a
mat33Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
2]
        m00 :: Holor a
m00 = Holor a
mat33Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2],[Int
1,Int
2]]
        m01 :: Holor a
m01 = Holor a
mat33Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2],[Int
0,Int
2]]
        m02 :: Holor a
m02 = Holor a
mat33Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2],[Int
0,Int
1]]

-- Laplace-Expansion by first row
-- | determinant of 4x4 matrix
det44 :: (Unbox a, Num a) => Holor a -> a
{-# INLINE det44 #-}
det44 :: Holor a -> a
det44 Holor a
mat44
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat44 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
4,Int
4] = a
h00a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 Holor a
m00) a -> a -> a
forall a. Num a => a -> a -> a
- a
h01a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 Holor a
m01) a -> a -> a
forall a. Num a => a -> a -> a
+ a
h02a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 Holor a
m02) a -> a -> a
forall a. Num a => a -> a -> a
- a
h03a -> a -> a
forall a. Num a => a -> a -> a
*(Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 Holor a
m03)
    | Bool
otherwise = [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char]
"det44 only possible for a holor with shape [4,4]")
    where
        h00 :: a
h00 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
0]
        h01 :: a
h01 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
1]
        h02 :: a
h02 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
2]
        h03 :: a
h03 = Holor a
mat44Holor a -> [Int] -> a
forall a. Unbox a => Holor a -> [Int] -> a
|!![Int
0,Int
3]
        m00 :: Holor a
m00 = Holor a
mat44Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2,Int
3],[Int
1,Int
2,Int
3]]
        m01 :: Holor a
m01 = Holor a
mat44Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2,Int
3],[Int
0,Int
2,Int
3]]
        m02 :: Holor a
m02 = Holor a
mat44Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2,Int
3],[Int
0,Int
1,Int
3]]
        m03 :: Holor a
m03 = Holor a
mat44Holor a -> [[Int]] -> Holor a
forall a. Unbox a => Holor a -> [[Int]] -> Holor a
|||!![[Int
1,Int
2,Int
3],[Int
0,Int
1,Int
2]]


{- TODO:
inv :: (Unbox a, Fractional a) => Holor a -> Holor a
{-# INLINE inv #-}
inv hlr = error "not implemented"
-}

-- Link: https://de.wikipedia.org/wiki/Inverse_Matrix#Explizite_Formeln
-- | inverse of 2x2 matrix
inv22 :: (Unbox a, Fractional a) => Holor a -> Holor a
{-# INLINE inv22 #-}
inv22 :: Holor a -> Holor a
inv22 Holor a
mat22
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat22 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
2,Int
2] = (a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det22 Holor a
mat22) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| (Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a
adjugate22 Holor a
mat22)
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char]
"inv22 only possible for a holor with shape [2,2]")

-- Link: https://de.wikipedia.org/wiki/Inverse_Matrix#Explizite_Formeln
-- | inverse of 3x3 matrix
inv33 :: (Unbox a, Fractional a) => Holor a -> Holor a
{-# INLINE inv33 #-}
inv33 :: Holor a -> Holor a
inv33 Holor a
mat33
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat33 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
3,Int
3] = (a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det33 Holor a
mat33) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| (Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a
adjugate33 Holor a
mat33)
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char]
"inv33 only possible for a holor with shape [3,3]")

-- Link: https://semath.info/src/inverse-cofactor-ex4.html
-- | inverse of 4x4 matrix
inv44 :: (Unbox a, Fractional a) => Holor a -> Holor a
{-# INLINE inv44 #-}
inv44 :: Holor a -> Holor a
inv44 Holor a
mat44
    | Holor a -> [Int]
forall a. Unbox a => Holor a -> [Int]
shape Holor a
mat44 [Int] -> [Int] -> Bool
forall a. Eq a => a -> a -> Bool
== [Int
4,Int
4] = (a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/Holor a -> a
forall a. (Unbox a, Num a) => Holor a -> a
det44 Holor a
mat44) a -> Holor a -> Holor a
forall a. (Unbox a, Num a) => a -> Holor a -> Holor a
*| (Holor a -> Holor a
forall a. (Unbox a, Num a) => Holor a -> Holor a
adjugate44 Holor a
mat44)
    | Bool
otherwise = [Char] -> Holor a
forall a. HasCallStack => [Char] -> a
error ([Char]
"inv44 only possible for a holor with shape [4,4]")


-- | select all indices of a holor with nonzero elements
selectNonzero :: (Unbox a, Num a, Eq a) => Holor a -> [[Int]]
{-# INLINE selectNonzero #-}
selectNonzero :: Holor a -> [[Int]]
selectNonzero Holor a
hlr = (a -> Bool) -> Holor a -> [[Int]]
forall a. Unbox a => (a -> Bool) -> Holor a -> [[Int]]
selectBy (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Holor a
hlr

-- | number of nonzero elements
countNonzero :: (Unbox a, Num a, Eq a) => Holor a -> Int
{-# INLINE countNonzero #-}
countNonzero :: Holor a -> Int
countNonzero Holor a
hlr = (a -> Bool) -> Holor a -> Int
forall a. Unbox a => (a -> Bool) -> Holor a -> Int
countBy (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/=a
0) Holor a
hlr