{-# LANGUAGE DuplicateRecordFields #-}
module Numeric.Integration.SimplexCubature
  (Result(..), Results(..), integrateOnSimplex, integrateOnSimplex')
  where
import           Data.Array.Unboxed                          (UArray, array)
import           Data.Array.Unsafe                           (unsafeThaw)
import qualified Data.Vector.Unboxed                         as UV
import           Numeric.Integration.Simplex.Simplex
import           Numeric.Integration.SimplexCubature.Internal

data Results = Results
  { Results -> [Double]
values         :: [Double]
  , Results -> [Double]
errorEstimates :: [Double]
  , Results -> Int
evaluations    :: Int
  , Results -> Bool
success        :: Bool
  } deriving Int -> Results -> ShowS
[Results] -> ShowS
Results -> String
(Int -> Results -> ShowS)
-> (Results -> String) -> ([Results] -> ShowS) -> Show Results
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Results] -> ShowS
$cshowList :: [Results] -> ShowS
show :: Results -> String
$cshow :: Results -> String
showsPrec :: Int -> Results -> ShowS
$cshowsPrec :: Int -> Results -> ShowS
Show

data Result = Result
  { Result -> Double
value         :: Double
  , Result -> Double
errorEstimate :: Double
  , Result -> Int
evaluations   :: Int
  , Result -> Bool
success       :: Bool
  } deriving Int -> Result -> ShowS
[Result] -> ShowS
Result -> String
(Int -> Result -> ShowS)
-> (Result -> String) -> ([Result] -> ShowS) -> Show Result
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Result] -> ShowS
$cshowList :: [Result] -> ShowS
show :: Result -> String
$cshow :: Result -> String
showsPrec :: Int -> Result -> ShowS
$cshowsPrec :: Int -> Result -> ShowS
Show

simplicesToArray :: Simplices -> IO IO3dArray
simplicesToArray :: Simplices -> IO IO3dArray
simplicesToArray Simplices
simplices = do
  let dim :: Int
dim = [Double] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[Double]] -> [Double]
forall a. [a] -> a
head (Simplices -> [[Double]]
forall a. [a] -> a
head Simplices
simplices))
      nsimplices :: Int
nsimplices = Simplices -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length Simplices
simplices
      assocList :: [((Int, Int, Int), Double)]
assocList = ([Int] -> ((Int, Int, Int), Double))
-> [[Int]] -> [((Int, Int, Int), Double)]
forall a b. (a -> b) -> [a] -> [b]
map (\[Int
i,Int
j,Int
k] -> ((Int
i,Int
j,Int
k), (Simplices
simplicesSimplices -> Int -> [[Double]]
forall a. [a] -> Int -> a
!!(Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))[[Double]] -> Int -> [Double]
forall a. [a] -> Int -> a
!!(Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)[Double] -> Int -> Double
forall a. [a] -> Int -> a
!!(Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)))
                      ([[Int]] -> [[Int]]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [[Int
1..Int
dim], [Int
1..(Int
dimInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)], [Int
1..Int
nsimplices]])
      arr :: UArray (Int, Int, Int) Double
arr = ((Int, Int, Int), (Int, Int, Int))
-> [((Int, Int, Int), Double)] -> UArray (Int, Int, Int) Double
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
(i, i) -> [(i, e)] -> a i e
array ((Int
1,Int
1,Int
1),(Int
dim,Int
dimInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1,Int
nsimplices)) [((Int, Int, Int), Double)]
assocList
            :: UArray (Int,Int,Int) Double
  UArray (Int, Int, Int) Double -> IO IO3dArray
forall i (a :: * -> * -> *) e (b :: * -> * -> *) (m :: * -> *).
(Ix i, IArray a e, MArray b e m) =>
a i e -> m (b i e)
unsafeThaw UArray (Int, Int, Int) Double
arr

-- | Integral of a vector-valued function over an union of simplices.

integrateOnSimplex
    :: (VectorD -> VectorD)   -- integrand

    -> Simplices              -- domain

    -> Int                    -- number of components

    -> Int                    -- maximum number of evaluations

    -> Double                 -- desired absolute error

    -> Double                 -- desired relative error

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

    -> IO Results             -- integral, error, evaluations, success

integrateOnSimplex :: (VectorD -> VectorD)
-> Simplices -> Int -> Int -> Double -> Double -> Int -> IO Results
integrateOnSimplex VectorD -> VectorD
f Simplices
s Int
ncomp Int
maxevals Double
absError Double
relError Int
rule = do
  let n :: Int
n = [[Double]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Simplices -> [[Double]]
forall a. [a] -> a
head Simplices
s) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
  if Simplices -> Bool
isValidSimplices Simplices
s
    then do
      IO3dArray
v <- Simplices -> IO IO3dArray
simplicesToArray Simplices
s
      (VectorD
vals, VectorD
errors, Int
nevals, Bool
fl) <-
        Int
-> Int
-> Int
-> (VectorD -> VectorD)
-> Double
-> Double
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
adsimp Int
n Int
ncomp Int
maxevals VectorD -> VectorD
f Double
absError Double
relError Int
rule IO3dArray
v
      Results -> IO Results
forall (m :: * -> *) a. Monad m => a -> m a
return (Results -> IO Results) -> Results -> IO Results
forall a b. (a -> b) -> a -> b
$ [Double] -> [Double] -> Int -> Bool -> Results
Results (VectorD -> [Double]
forall a. Unbox a => Vector a -> [a]
UV.toList VectorD
vals) (VectorD -> [Double]
forall a. Unbox a => Vector a -> [a]
UV.toList VectorD
errors) Int
nevals (Bool -> Bool
not Bool
fl)
    else String -> IO Results
forall a. HasCallStack => String -> a
error String
"invalid simplices"

-- | Integral of a real-valued function over an union of simplices.

integrateOnSimplex'
    :: (VectorD -> Double)    -- integrand

    -> Simplices              -- domain

    -> Int                    -- maximum number of evaluations

    -> Double                 -- desired absolute error

    -> Double                 -- desired relative error

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

    -> IO Result              -- integral, error, evaluations, success

integrateOnSimplex' :: (VectorD -> Double)
-> Simplices -> Int -> Double -> Double -> Int -> IO Result
integrateOnSimplex' VectorD -> Double
f Simplices
s Int
maxevals Double
absError Double
relError Int
rule = do
  let n :: Int
n = [[Double]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length (Simplices -> [[Double]]
forall a. [a] -> a
head Simplices
s) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
  if Simplices -> Bool
isValidSimplices Simplices
s
    then do
      IO3dArray
v <- Simplices -> IO IO3dArray
simplicesToArray Simplices
s
      (VectorD
val, VectorD
err, Int
nevals, Bool
fl) <-
        Int
-> Int
-> Int
-> (VectorD -> VectorD)
-> Double
-> Double
-> Int
-> IO3dArray
-> IO (VectorD, VectorD, Int, Bool)
adsimp Int
n Int
1 Int
maxevals (Double -> VectorD
forall a. Unbox a => a -> Vector a
UV.singleton (Double -> VectorD) -> (VectorD -> Double) -> VectorD -> VectorD
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VectorD -> Double
f) Double
absError Double
relError Int
rule IO3dArray
v
      Result -> IO Result
forall (m :: * -> *) a. Monad m => a -> m a
return (Result -> IO Result) -> Result -> IO Result
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Int -> Bool -> Result
Result (VectorD -> Double
forall a. Unbox a => Vector a -> a
UV.head VectorD
val) (VectorD -> Double
forall a. Unbox a => Vector a -> a
UV.head VectorD
err) Int
nevals (Bool -> Bool
not Bool
fl)
    else String -> IO Result
forall a. HasCallStack => String -> a
error String
"invalid simplices"