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
dim Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Bool -> Bool -> Bool
&&
    ([Double] -> Bool) -> Simplex -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
dim) (Int -> Bool) -> ([Double] -> Int) -> [Double] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [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
&&
    (Simplex -> Bool) -> Simplices -> 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)) (Int -> Bool) -> (Simplex -> Int) -> Simplex -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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 -- rq: tu calcules le fact à chaque fois
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 -- not used
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)