module Numeric.Integration.SphericalSimplexCubature.Internal
(orthants, SphericalSimplex, transformedIntegrand)
where
import Data.List.Index (imap)
import Data.Matrix (detLU, diagonalList, fromLists,
minorMatrix, toLists, zero, (<->))
import Data.Vector.Unboxed (Vector)
import qualified Data.Vector.Unboxed as V
type SphericalSimplex = [[Double]]
orthants :: Int -> [SphericalSimplex]
orthants :: Int -> [SphericalSimplex]
orthants Int
n = forall a. [a] -> [a]
reverse forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Matrix a -> [[a]]
toLists forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Int -> a -> [a] -> Matrix a
diagonalList Int
n Double
0) (forall {t} {a}. (Eq t, Num t, Num a) => t -> [[a]]
pm Int
n)
where pm :: t -> [[a]]
pm t
2 = [[a
i, a
j] | a
i <- [-a
1, a
1], a
j <- [-a
1, a
1]]
pm t
k = [a
i forall a. a -> [a] -> [a]
: [a]
l | a
i <- [-a
1, a
1], [a]
l <- t -> [[a]]
pm (t
kforall a. Num a => a -> a -> a
-t
1)]
norm2 :: [Double] -> Double
norm2 :: [Double] -> Double
norm2 [Double]
v = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [Double]
v [Double]
v
dotproduct :: [Double] -> [Double] -> Double
dotproduct :: [Double] -> [Double] -> Double
dotproduct [Double]
a [Double]
b = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) [Double]
a [Double]
b
scalarTimesList :: Double -> [Double] -> [Double]
scalarTimesList :: Double -> [Double] -> [Double]
scalarTimesList Double
lambda = forall a b. (a -> b) -> [a] -> [b]
map (forall a. Num a => a -> a -> a
* Double
lambda)
f :: SphericalSimplex -> [Double] -> [Double]
f :: SphericalSimplex -> [Double] -> [Double]
f SphericalSimplex
vertices [Double]
stu = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(+)) (SphericalSimplex
verticesforall a. [a] -> Int -> a
!!Int
0) SphericalSimplex
terms
where
w :: SphericalSimplex
w = forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
subtract (SphericalSimplex
verticesforall a. [a] -> Int -> a
!!Int
0)) (forall a. [a] -> [a]
tail SphericalSimplex
vertices)
terms :: SphericalSimplex
terms = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> [Double] -> [Double]
scalarTimesList [Double]
stu SphericalSimplex
w
g :: SphericalSimplex -> [Double] -> [Double]
g :: SphericalSimplex -> [Double] -> [Double]
g SphericalSimplex
vertices [Double]
stu = Double -> [Double] -> [Double]
scalarTimesList (Double
1 forall a. Fractional a => a -> a -> a
/ [Double] -> Double
norm2 [Double]
fstu) [Double]
fstu
where fstu :: [Double]
fstu = SphericalSimplex -> [Double] -> [Double]
f SphericalSimplex
vertices [Double]
stu
dg :: SphericalSimplex -> [Double] -> [[Double]]
dg :: SphericalSimplex -> [Double] -> SphericalSimplex
dg SphericalSimplex
vertices [Double]
stu = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
subtract) SphericalSimplex
fviv1 SphericalSimplex
nviv1
where
fstu :: [Double]
fstu = SphericalSimplex -> [Double] -> [Double]
f SphericalSimplex
vertices [Double]
stu
invn :: Double
invn = Double
1 forall a. Fractional a => a -> a -> a
/ [Double] -> Double
norm2 [Double]
fstu
invn3 :: Double
invn3 = Double
invnforall a. Num a => a -> a -> a
*Double
invnforall a. Num a => a -> a -> a
*Double
invn
viv1 :: SphericalSimplex
viv1 = forall a b. (a -> b) -> [a] -> [b]
map (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
subtract (forall a. [a] -> a
head SphericalSimplex
vertices)) (forall a. [a] -> [a]
tail SphericalSimplex
vertices)
nviv1 :: SphericalSimplex
nviv1 = forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scalarTimesList Double
invn) SphericalSimplex
viv1
dpi :: [Double]
dpi = forall a b. (a -> b) -> [a] -> [b]
map ((forall a. Num a => a -> a -> a
*Double
invn3) forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double] -> Double
dotproduct [Double]
fstu) SphericalSimplex
viv1
fviv1 :: SphericalSimplex
fviv1 = forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
`scalarTimesList` [Double]
fstu) [Double]
dpi
extProduct :: [[Double]] -> [Double]
extProduct :: SphericalSimplex -> [Double]
extProduct SphericalSimplex
vectors =
forall a b. (Int -> a -> b) -> [a] -> [b]
imap (\Int
i Matrix Double
mat -> (if forall a. Integral a => a -> Bool
even Int
i then Double
1 else -Double
1) forall a. Num a => a -> a -> a
* forall a. (Ord a, Fractional a) => Matrix a -> a
detLU Matrix Double
mat) [Matrix Double]
minorMatrices
where
dim :: Int
dim = forall (t :: * -> *) a. Foldable t => t a -> Int
length SphericalSimplex
vectors forall a. Num a => a -> a -> a
+ Int
1
matrix :: Matrix Double
matrix = forall a. Num a => Int -> Int -> Matrix a
zero Int
1 Int
dim forall a. Matrix a -> Matrix a -> Matrix a
<-> forall a. [[a]] -> Matrix a
fromLists SphericalSimplex
vectors
minorMatrices :: [Matrix Double]
minorMatrices = forall a b. (a -> b) -> [a] -> [b]
map (\Int
j -> forall a. Int -> Int -> Matrix a -> Matrix a
minorMatrix Int
1 Int
j Matrix Double
matrix) [Int
1 .. Int
dim]
sigma :: SphericalSimplex -> [Double] -> Double
sigma :: SphericalSimplex -> [Double] -> Double
sigma SphericalSimplex
ssimplex [Double]
stu = [Double] -> Double
norm2 forall a b. (a -> b) -> a -> b
$ SphericalSimplex -> [Double]
extProduct (SphericalSimplex -> [Double] -> SphericalSimplex
dg SphericalSimplex
ssimplex [Double]
stu)
transformedIntegrand :: SphericalSimplex -> ([Double] -> Double)
-> (Vector Double -> Double)
transformedIntegrand :: SphericalSimplex -> ([Double] -> Double) -> Vector Double -> Double
transformedIntegrand SphericalSimplex
ssimplex [Double] -> Double
integrand Vector Double
stu =
let stul :: [Double]
stul = forall a. Unbox a => Vector a -> [a]
V.toList Vector Double
stu in
SphericalSimplex -> [Double] -> Double
sigma SphericalSimplex
ssimplex [Double]
stul forall a. Num a => a -> a -> a
* [Double] -> Double
integrand (SphericalSimplex -> [Double] -> [Double]
g SphericalSimplex
ssimplex [Double]
stul)