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]] -- square [v1, v2, v3, v4]

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)