{-# 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 ( isValidSimplices, Simplices )
import           Numeric.Integration.SimplexCubature.Internal  ( VectorD, IO3dArray, adsimp )

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
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
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 = forall (t :: * -> *) a. Foldable t => t a -> Int
length (forall a. [a] -> a
head (forall a. [a] -> a
head Simplices
simplices))
      nsimplices :: Int
nsimplices = forall (t :: * -> *) a. Foldable t => t a -> Int
length Simplices
simplices
      assocList :: [((Int, Int, Int), Double)]
assocList = forall a b. (a -> b) -> [a] -> [b]
map (\[Int
i, Int
j, Int
k] -> ((Int
i, Int
j, Int
k), (Simplices
simplicesforall a. [a] -> Int -> a
!!(Int
kforall a. Num a => a -> a -> a
-Int
1))forall a. [a] -> Int -> a
!!(Int
jforall a. Num a => a -> a -> a
-Int
1)forall a. [a] -> Int -> a
!!(Int
iforall a. Num a => a -> a -> a
-Int
1)))
                      (forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [[Int
1 .. Int
dim], [Int
1 .. (Int
dimforall a. Num a => a -> a -> a
+Int
1)], [Int
1 .. Int
nsimplices]])
      arr :: UArray (Int, Int, Int) Double
arr = 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
dimforall a. Num a => a -> a -> a
+Int
1, Int
nsimplices)) [((Int, Int, Int), Double)]
assocList
            :: UArray (Int, Int, Int) Double
  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 of integration
    -> 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 = forall (t :: * -> *) a. Foldable t => t a -> Int
length (forall a. [a] -> a
head Simplices
s) 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
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ [Double] -> [Double] -> Int -> Bool -> Results
Results (forall a. Unbox a => Vector a -> [a]
UV.toList VectorD
vals) (forall a. Unbox a => Vector a -> [a]
UV.toList VectorD
errors) Int
nevals (Bool -> Bool
not Bool
fl)
    else 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 of integration
    -> 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 = forall (t :: * -> *) a. Foldable t => t a -> Int
length (forall a. [a] -> a
head Simplices
s) 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 (forall a. Unbox a => a -> Vector a
UV.singleton forall b c a. (b -> c) -> (a -> b) -> a -> c
. VectorD -> Double
f) Double
absError Double
relError Int
rule IO3dArray
v
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Double -> Double -> Int -> Bool -> Result
Result (forall a. Unbox a => Vector a -> a
UV.head VectorD
val) (forall a. Unbox a => Vector a -> a
UV.head VectorD
err) Int
nevals (Bool -> Bool
not Bool
fl)
    else forall a. HasCallStack => String -> a
error String
"invalid simplices"