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 = [SphericalSimplex] -> [SphericalSimplex]
forall a. [a] -> [a]
reverse ([SphericalSimplex] -> [SphericalSimplex])
-> [SphericalSimplex] -> [SphericalSimplex]
forall a b. (a -> b) -> a -> b
$ ([Double] -> SphericalSimplex)
-> SphericalSimplex -> [SphericalSimplex]
forall a b. (a -> b) -> [a] -> [b]
map (Matrix Double -> SphericalSimplex
forall a. Matrix a -> [[a]]
toLists (Matrix Double -> SphericalSimplex)
-> ([Double] -> Matrix Double) -> [Double] -> SphericalSimplex
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double -> [Double] -> Matrix Double
forall a. Int -> a -> [a] -> Matrix a
diagonalList Int
n Double
0) (Int -> SphericalSimplex
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 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a]
l | a
i <- [-a
1,a
1], [a]
l <- t -> [[a]]
pm (t
kt -> t -> t
forall a. Num a => a -> a -> a
-t
1)]

norm2 :: [Double] -> Double
norm2 :: [Double] -> Double
norm2 [Double]
v = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ [Double] -> Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Double] -> Double) -> [Double] -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
v [Double]
v

dotproduct :: [Double] -> [Double] -> Double
dotproduct :: [Double] -> [Double] -> Double
dotproduct [Double]
a [Double]
b = [Double] -> Double
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Double] -> Double) -> [Double] -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
a [Double]
b

scalarTimesList :: Double -> [Double] -> [Double]
scalarTimesList :: Double -> [Double] -> [Double]
scalarTimesList Double
lambda = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lambda)

f :: SphericalSimplex -> [Double] -> [Double]
f :: SphericalSimplex -> [Double] -> [Double]
f SphericalSimplex
vertices [Double]
stu = ([Double] -> [Double] -> [Double])
-> [Double] -> SphericalSimplex -> [Double]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)) (SphericalSimplex
verticesSphericalSimplex -> Int -> [Double]
forall a. [a] -> Int -> a
!!Int
0) SphericalSimplex
terms
  where
    w :: SphericalSimplex
w = ([Double] -> [Double]) -> SphericalSimplex -> SphericalSimplex
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract (SphericalSimplex
verticesSphericalSimplex -> Int -> [Double]
forall a. [a] -> Int -> a
!!Int
0)) (SphericalSimplex -> SphericalSimplex
forall a. [a] -> [a]
tail SphericalSimplex
vertices)
    terms :: SphericalSimplex
terms = (Double -> [Double] -> [Double])
-> [Double] -> SphericalSimplex -> SphericalSimplex
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 Double -> Double -> Double
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 = ([Double] -> [Double] -> [Double])
-> SphericalSimplex -> SphericalSimplex -> SphericalSimplex
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
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 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ [Double] -> Double
norm2 [Double]
fstu
    invn3 :: Double
invn3 = Double
invnDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invnDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invn
    viv1 :: SphericalSimplex
viv1 = ([Double] -> [Double]) -> SphericalSimplex -> SphericalSimplex
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract (SphericalSimplex -> [Double]
forall a. [a] -> a
head SphericalSimplex
vertices)) (SphericalSimplex -> SphericalSimplex
forall a. [a] -> [a]
tail SphericalSimplex
vertices)
    nviv1 :: SphericalSimplex
nviv1 = ([Double] -> [Double]) -> SphericalSimplex -> SphericalSimplex
forall a b. (a -> b) -> [a] -> [b]
map (Double -> [Double] -> [Double]
scalarTimesList Double
invn) SphericalSimplex
viv1
    dpi :: [Double]
dpi = ([Double] -> Double) -> SphericalSimplex -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
invn3) (Double -> Double) -> ([Double] -> Double) -> [Double] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double] -> Double
dotproduct [Double]
fstu) SphericalSimplex
viv1
    fviv1 :: SphericalSimplex
fviv1 = (Double -> [Double]) -> [Double] -> SphericalSimplex
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 =
  (Int -> Matrix Double -> Double) -> [Matrix Double] -> [Double]
forall a b. (Int -> a -> b) -> [a] -> [b]
imap (\Int
i Matrix Double
mat -> (if Int -> Bool
forall a. Integral a => a -> Bool
even Int
i then Double
1 else -Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Matrix Double -> Double
forall a. (Ord a, Fractional a) => Matrix a -> a
detLU Matrix Double
mat) [Matrix Double]
minorMatrices
  where
    dim :: Int
dim = SphericalSimplex -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length SphericalSimplex
vectors Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
    matrix :: Matrix Double
matrix = Int -> Int -> Matrix Double
forall a. Num a => Int -> Int -> Matrix a
zero Int
1 Int
dim Matrix Double -> Matrix Double -> Matrix Double
forall a. Matrix a -> Matrix a -> Matrix a
<-> SphericalSimplex -> Matrix Double
forall a. [[a]] -> Matrix a
fromLists SphericalSimplex
vectors
    minorMatrices :: [Matrix Double]
minorMatrices = (Int -> Matrix Double) -> [Int] -> [Matrix Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Int
j -> Int -> Int -> Matrix Double -> Matrix Double
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 ([Double] -> Double) -> [Double] -> Double
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 = Vector Double -> [Double]
forall a. Unbox a => Vector a -> [a]
V.toList Vector Double
stu in
  SphericalSimplex -> [Double] -> Double
sigma SphericalSimplex
ssimplex [Double]
stul Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
integrand (SphericalSimplex -> [Double] -> [Double]
g SphericalSimplex
ssimplex [Double]
stul)