module Numeric.Integration.IntegratePolynomialOnSimplex where
import           Algebra.Ring                   ( C )
import           Data.List                      ( transpose )
import           Data.Matrix                    ( detLU
                                                , fromLists
                                                )
import           Math.Algebra.Hspray            ( (*^)
                                                , Spray
                                                , (^+^)
                                                , bombieriSpray
                                                , composeSpray
                                                , constantSpray
                                                , lone
                                                , toList
                                                )

-- | Exact integral of a polynomial over a simplex
integratePolynomialOnSimplex
  :: (C a, Fractional a, Ord a) 
  => Spray a -- ^ polynomial to be integrated
  -> [[a]]   -- ^ simplex to integrate over
  -> a
integratePolynomialOnSimplex :: forall a. (C a, Fractional a, Ord a) => Spray a -> [[a]] -> a
integratePolynomialOnSimplex Spray a
p [[a]]
simplex =
  a
s forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
abs (forall a. (Ord a, Fractional a) => Matrix a -> a
detLU forall a b. (a -> b) -> a -> b
$ forall a. [[a]] -> Matrix a
fromLists [[a]]
b) forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int
2 .. Int
n])
 where
  v :: [a]
v            = forall a. [a] -> a
last [[a]]
simplex
  n :: Int
n            = forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
v
  b :: [[a]]
b            = forall a b. (a -> b) -> [a] -> [b]
map (\[a]
column -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [a]
column [a]
v) (forall a. Int -> [a] -> [a]
take Int
n [[a]]
simplex)
  vb :: [(a, [a])]
vb           = forall a b. [a] -> [b] -> [(a, b)]
zip [a]
v (forall a. [[a]] -> [[a]]
transpose [[a]]
b)
  variables :: [Spray a]
variables    = forall a b. (a -> b) -> [a] -> [b]
map forall a. C a => Int -> Spray a
lone [Int
1 .. Int
n]
  newvariables :: [Spray a]
newvariables = forall a b. (a -> b) -> [a] -> [b]
map
    (\(a
vi, [a]
bi) ->
      (forall a. (C a, Eq a) => a -> Spray a
constantSpray a
vi) forall a. (C a, Eq a) => Spray a -> Spray a -> Spray a
^+^ forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall a. (C a, Eq a) => Spray a -> Spray a -> Spray a
(^+^) (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. (C a, Eq a) => a -> Spray a -> Spray a
(*^) [a]
bi [Spray a]
variables)
    )
    [(a, [a])]
vb
  q :: Spray a
q      = forall a. (C a, Eq a) => Spray a -> [Spray a] -> Spray a
composeSpray Spray a
p [Spray a]
newvariables
  qterms :: [([Int], a)]
qterms = forall a. Spray a -> [([Int], a)]
toList forall a b. (a -> b) -> a -> b
$ forall a. C a => Spray a -> Spray a
bombieriSpray Spray a
q
  s :: a
s      = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall {b} {t :: * -> *}.
(Fractional b, Foldable t) =>
(t Int, b) -> b
f [([Int], a)]
qterms
   where
    f :: (t Int, b) -> b
f (t Int
exponents, b
coef) = if Int
d forall a. Eq a => a -> a -> Bool
== Int
0
      then b
coef
      else b
coef forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int
n forall a. Num a => a -> a -> a
+ Int
1 .. Int
n forall a. Num a => a -> a -> a
+ Int
d])
      where d :: Int
d = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum t Int
exponents