{-|
Module      : Numeric.Integration.PolyhedralCubature
Description : Multiple integration over convex polytopes.
Copyright   : (c) Stéphane Laurent, 2023-2024
License     : GPL-3
Maintainer  : laurent_step@outlook.fr

Evaluation of integrals over a convex polytope. See README for examples.
-}
module Numeric.Integration.PolyhedralCubature
  ( 
    Result(..)
  , Results(..)
  , Constraint(..)
  , VectorD
  , integrateOnPolytopeN
  , integrateOnPolytope
  , integrateOnPolytopeN'
  , integrateOnPolytope'
  , integratePolynomialOnPolytope
  , integratePolynomialOnPolytope'
  )
  where
import Algebra.Ring                                     ( C )
import qualified Data.IntMap.Strict                     as IM
import Data.Vector.Unboxed                              ( Vector )
import Geometry.Delaunay                                ( delaunay
                                                        , getDelaunayTiles 
                                                        )
import Geometry.VertexEnum                              ( Constraint(..)
                                                        , vertexenum 
                                                        )
import Math.Algebra.Hspray                              ( Spray )
import Numeric.Integration.IntegratePolynomialOnSimplex ( integratePolynomialOnSimplex )
import Numeric.Integration.SimplexCubature              ( Result(..)
                                                        , Results(..)
                                                        , integrateOnSimplex
                                                        , integrateOnSimplex'
                                                        )

type VectorD = Vector Double

-- | Integral of a multivariate function over a convex polytope given by its vertices.

integrateOnPolytopeN
    :: (VectorD -> VectorD)   -- ^ integrand

    -> [[Double]]             -- ^ vertices of the polytope

    -> Int                    -- ^ number of components of the integrand

    -> Int                    -- ^ maximum number of evaluations

    -> Double                 -- ^ desired absolute error

    -> Double                 -- ^ desired relative error

    -> Int                    -- ^ integration rule: 1, 2, 3 or 4

    -> IO Results             -- ^ values, error estimate, evaluations, success

integrateOnPolytopeN :: (VectorD -> VectorD)
-> [[Double]]
-> Int
-> Int
-> Double
-> Double
-> Int
-> IO Results
integrateOnPolytopeN VectorD -> VectorD
f [[Double]]
vertices Int
dim Int
maxevals Double
abserr Double
relerr Int
rule = do 
  Tessellation
tessellation <- [[Double]] -> Bool -> Bool -> Maybe Double -> IO Tessellation
delaunay [[Double]]
vertices Bool
False Bool
False Maybe Double
forall a. Maybe a
Nothing
  let simplices :: [[[Double]]]
simplices = (IntMap [Double] -> [[Double]])
-> [IntMap [Double]] -> [[[Double]]]
forall a b. (a -> b) -> [a] -> [b]
map IntMap [Double] -> [[Double]]
forall a. IntMap a -> [a]
IM.elems (Tessellation -> [IntMap [Double]]
getDelaunayTiles Tessellation
tessellation)
  (VectorD -> VectorD)
-> [[[Double]]]
-> Int
-> Int
-> Double
-> Double
-> Int
-> IO Results
integrateOnSimplex VectorD -> VectorD
f [[[Double]]]
simplices Int
dim Int
maxevals Double
abserr Double
relerr Int
rule

-- | Integral of a real-valued function over a convex polytope given by its vertices.

integrateOnPolytope
    :: (VectorD -> Double)    -- ^ integrand

    -> [[Double]]             -- ^ vertices of the polytope

    -> Int                    -- ^ maximum number of evaluations

    -> Double                 -- ^ desired absolute error

    -> Double                 -- ^ desired relative error

    -> Int                    -- ^ integration rule: 1, 2, 3 or 4

    -> IO Result              -- ^ values, error estimate, evaluations, success

integrateOnPolytope :: (VectorD -> Double)
-> [[Double]] -> Int -> Double -> Double -> Int -> IO Result
integrateOnPolytope VectorD -> Double
f [[Double]]
vertices Int
maxevals Double
abserr Double
relerr Int
rule = do 
  Tessellation
tessellation <- [[Double]] -> Bool -> Bool -> Maybe Double -> IO Tessellation
delaunay [[Double]]
vertices Bool
True Bool
False Maybe Double
forall a. Maybe a
Nothing
  let simplices :: [[[Double]]]
simplices = (IntMap [Double] -> [[Double]])
-> [IntMap [Double]] -> [[[Double]]]
forall a b. (a -> b) -> [a] -> [b]
map IntMap [Double] -> [[Double]]
forall a. IntMap a -> [a]
IM.elems (Tessellation -> [IntMap [Double]]
getDelaunayTiles Tessellation
tessellation)
  (VectorD -> Double)
-> [[[Double]]] -> Int -> Double -> Double -> Int -> IO Result
integrateOnSimplex' VectorD -> Double
f [[[Double]]]
simplices Int
maxevals Double
abserr Double
relerr Int
rule

-- | Integral of a multivariate function over a convex polytope given by linear inequalities.

integrateOnPolytopeN'
    :: Real a
    => (VectorD -> VectorD)   -- ^ integrand

    -> [Constraint a]         -- ^ linear inequalities defining the polytope

    -> Int                    -- ^ number of components of the integrand

    -> Int                    -- ^ maximum number of evaluations

    -> Double                 -- ^ desired absolute error

    -> Double                 -- ^ desired relative error

    -> Int                    -- ^ integration rule: 1, 2, 3 or 4

    -> IO Results             -- ^ values, error estimate, evaluations, success

integrateOnPolytopeN' :: forall a.
Real a =>
(VectorD -> VectorD)
-> [Constraint a]
-> Int
-> Int
-> Double
-> Double
-> Int
-> IO Results
integrateOnPolytopeN' VectorD -> VectorD
f [Constraint a]
constraints Int
dim Int
maxevals Double
abserr Double
relerr Int
rule = do 
  [[Double]]
vertices <- [Constraint a] -> Maybe [Double] -> IO [[Double]]
forall a.
Real a =>
[Constraint a] -> Maybe [Double] -> IO [[Double]]
vertexenum [Constraint a]
constraints Maybe [Double]
forall a. Maybe a
Nothing
  (VectorD -> VectorD)
-> [[Double]]
-> Int
-> Int
-> Double
-> Double
-> Int
-> IO Results
integrateOnPolytopeN VectorD -> VectorD
f [[Double]]
vertices Int
dim Int
maxevals Double
abserr Double
relerr Int
rule

-- | Integral of a scalar-valued function over a convex polytope given by linear inequalities.

integrateOnPolytope'
    :: Real a
    => (VectorD -> Double)    -- ^ integrand

    -> [Constraint a]         -- ^ linear inequalities defining the polytope

    -> Int                    -- ^ maximum number of evaluations

    -> Double                 -- ^ desired absolute error

    -> Double                 -- ^ desired relative error

    -> Int                    -- ^ integration rule: 1, 2, 3 or 4

    -> IO Result              -- ^ values, error estimate, evaluations, success

integrateOnPolytope' :: forall a.
Real a =>
(VectorD -> Double)
-> [Constraint a] -> Int -> Double -> Double -> Int -> IO Result
integrateOnPolytope' VectorD -> Double
f [Constraint a]
constraints Int
maxevals Double
abserr Double
relerr Int
rule = do 
  [[Double]]
vertices <- [Constraint a] -> Maybe [Double] -> IO [[Double]]
forall a.
Real a =>
[Constraint a] -> Maybe [Double] -> IO [[Double]]
vertexenum [Constraint a]
constraints Maybe [Double]
forall a. Maybe a
Nothing
  (VectorD -> Double)
-> [[Double]] -> Int -> Double -> Double -> Int -> IO Result
integrateOnPolytope VectorD -> Double
f [[Double]]
vertices Int
maxevals Double
abserr Double
relerr Int
rule

delaunay' :: Real a => [[a]] -> IO [[[a]]]
delaunay' :: forall a. Real a => [[a]] -> IO [[[a]]]
delaunay' [[a]]
points = do
  let points' :: [[Double]]
points' = ([a] -> [Double]) -> [[a]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((a -> Double) -> [a] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map a -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac) [[a]]
points
  Tessellation
tessellation <- [[Double]] -> Bool -> Bool -> Maybe Double -> IO Tessellation
delaunay [[Double]]
points' Bool
True Bool
False Maybe Double
forall a. Maybe a
Nothing
  let indices :: [[Int]]
indices = (IntMap [Double] -> [Int]) -> [IntMap [Double]] -> [[Int]]
forall a b. (a -> b) -> [a] -> [b]
map IntMap [Double] -> [Int]
forall a. IntMap a -> [Int]
IM.keys (Tessellation -> [IntMap [Double]]
getDelaunayTiles Tessellation
tessellation)
  [[[a]]] -> IO [[[a]]]
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ([[[a]]] -> IO [[[a]]]) -> [[[a]]] -> IO [[[a]]]
forall a b. (a -> b) -> a -> b
$ ([Int] -> [[a]]) -> [[Int]] -> [[[a]]]
forall a b. (a -> b) -> [a] -> [b]
map ((Int -> [a]) -> [Int] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map ([[a]]
points [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!!)) [[Int]]
indices

-- | Integral of a polynomial over a convex polytope given by its vertices.

integratePolynomialOnPolytope
  :: (RealFrac a, C a)
  => Spray a -- ^ polynomial to be integrated

  -> [[a]]   -- ^ vertices of the polytope to integrate over

  -> IO a
integratePolynomialOnPolytope :: forall a. (RealFrac a, C a) => Spray a -> [[a]] -> IO a
integratePolynomialOnPolytope Spray a
p [[a]]
vertices = do
  [[[a]]]
simplices <- [[a]] -> IO [[[a]]]
forall a. Real a => [[a]] -> IO [[[a]]]
delaunay' [[a]]
vertices
  let integrals :: [a]
integrals = ([[a]] -> a) -> [[[a]]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (Spray a -> [[a]] -> a
forall a. (C a, Fractional a, Ord a) => Spray a -> [[a]] -> a
integratePolynomialOnSimplex Spray a
p) [[[a]]]
simplices 
  a -> IO a
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> a -> IO a
forall a b. (a -> b) -> a -> b
$ [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a]
integrals

-- | Integral of a polynomial over a convex polytope given by linear inequalities.

integratePolynomialOnPolytope'
  :: Spray Double        -- ^ polynomial to be integrated

  -> [Constraint Double] -- ^ linear inequalities defining the polytope

  -> IO Double
integratePolynomialOnPolytope' :: Spray Double -> [Constraint Double] -> IO Double
integratePolynomialOnPolytope' Spray Double
p [Constraint Double]
constraints = do
  [[Double]]
vertices <- [Constraint Double] -> Maybe [Double] -> IO [[Double]]
forall a.
Real a =>
[Constraint a] -> Maybe [Double] -> IO [[Double]]
vertexenum [Constraint Double]
constraints Maybe [Double]
forall a. Maybe a
Nothing
  Tessellation
tessellation <- [[Double]] -> Bool -> Bool -> Maybe Double -> IO Tessellation
delaunay [[Double]]
vertices Bool
True Bool
False Maybe Double
forall a. Maybe a
Nothing
  let simplices :: [[[Double]]]
simplices = (IntMap [Double] -> [[Double]])
-> [IntMap [Double]] -> [[[Double]]]
forall a b. (a -> b) -> [a] -> [b]
map IntMap [Double] -> [[Double]]
forall a. IntMap a -> [a]
IM.elems (Tessellation -> [IntMap [Double]]
getDelaunayTiles Tessellation
tessellation)
      integrals :: [Double]
integrals = ([[Double]] -> Double) -> [[[Double]]] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Spray Double -> [[Double]] -> Double
forall a. (C a, Fractional a, Ord a) => Spray a -> [[a]] -> a
integratePolynomialOnSimplex Spray Double
p) [[[Double]]]
simplices 
  Double -> IO Double
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Double -> IO Double) -> Double -> IO Double
forall a b. (a -> b) -> a -> b
$ [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Double]
integrals