module Numeric.Integration.Simplex.Simplex
where
import Data.Matrix (detLU, elementwiseUnsafe, fromLists)
type Simplex = [[Double]]
type Simplices = [Simplex]
isValidSimplex :: Simplex -> Bool
isValidSimplex :: Simplex -> Bool
isValidSimplex Simplex
simplex =
(Simplex -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Simplex
simplex Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
dimInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) Bool -> Bool -> Bool
&&
(Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
dim) (([Double] -> Int) -> Simplex -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map [Double] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Simplex -> Simplex
forall a. [a] -> [a]
tail Simplex
simplex))
where dim :: Int
dim = [Double] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Simplex -> [Double]
forall a. [a] -> a
head Simplex
simplex)
isValidSimplices :: Simplices -> Bool
isValidSimplices :: Simplices -> Bool
isValidSimplices Simplices
simplices =
(Simplex -> Bool) -> Simplices -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Simplex -> Bool
isValidSimplex Simplices
simplices Bool -> Bool -> Bool
&&
(Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Simplex -> Int
forall (t :: * -> *) a. Foldable t => [t a] -> Int
spaceDim (Simplices -> Simplex
forall a. [a] -> a
head Simplices
simplices)) ((Simplex -> Int) -> Simplices -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Simplex -> Int
forall (t :: * -> *) a. Foldable t => [t a] -> Int
spaceDim (Simplices -> Simplices
forall a. [a] -> [a]
tail Simplices
simplices))
where spaceDim :: [t a] -> Int
spaceDim [t a]
simplex = t a -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([t a] -> t a
forall a. [a] -> a
head [t a]
simplex)
canonicalSimplex :: Int -> Simplex
canonicalSimplex :: Int -> Simplex
canonicalSimplex Int
dim =
Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate Int
dim Double
0 [Double] -> Simplex -> Simplex
forall a. a -> [a] -> [a]
:
(Int -> [Double]) -> [Int] -> Simplex
forall a b. (a -> b) -> [a] -> [b]
map (\Int
v -> (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral(Int -> Double) -> (Int -> Int) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.Bool -> Int
forall a. Enum a => a -> Int
fromEnum(Bool -> Int) -> (Int -> Bool) -> Int -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
.(Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
v)) [Int
1..Int
dim]) [Int
1..Int
dim]
simplexVolume :: Simplex -> Double
simplexVolume :: Simplex -> Double
simplexVolume Simplex
s = Double -> Double
forall a. Num a => a -> a
abs (Matrix Double -> Double
forall a. (Ord a, Fractional a) => Matrix a -> a
detLU Matrix Double
v) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int] -> Int
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int
1..Int
n])
where n :: Int
n = Simplex -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Simplex
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
m1 :: Matrix Double
m1 = Simplex -> Matrix Double
forall a. [[a]] -> Matrix a
fromLists (Simplex -> Simplex
forall a. [a] -> [a]
tail Simplex
s)
m2 :: Matrix Double
m2 = Simplex -> Matrix Double
forall a. [[a]] -> Matrix a
fromLists (Simplex -> Matrix Double) -> Simplex -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Int -> [Double] -> Simplex
forall a. Int -> a -> [a]
replicate Int
n (Simplex -> [Double]
forall a. [a] -> a
head Simplex
s)
v :: Matrix Double
v = (Double -> Double -> Double)
-> Matrix Double -> Matrix Double -> Matrix Double
forall a b c. (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
elementwiseUnsafe (-) Matrix Double
m1 Matrix Double
m2
jacobian :: Simplex -> Double
jacobian :: Simplex -> Double
jacobian Simplex
s = Double -> Double
forall a. Num a => a -> a
abs (Matrix Double -> Double
forall a. (Ord a, Fractional a) => Matrix a -> a
detLU ((Double -> Double -> Double)
-> Matrix Double -> Matrix Double -> Matrix Double
forall a b c. (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
elementwiseUnsafe (-) Matrix Double
m1 Matrix Double
m2))
where m1 :: Matrix Double
m1 = Simplex -> Matrix Double
forall a. [[a]] -> Matrix a
fromLists (Simplex -> Simplex
forall a. [a] -> [a]
tail Simplex
s)
m2 :: Matrix Double
m2 = Simplex -> Matrix Double
forall a. [[a]] -> Matrix a
fromLists (Simplex -> Matrix Double) -> Simplex -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Int -> [Double] -> Simplex
forall a. Int -> a -> [a]
replicate (Simplex -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Simplex
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Simplex -> [Double]
forall a. [a] -> a
head Simplex
s)