{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
module Data.Massiv.Array.Numeric.Integral
(
midpointRule
, midpointStencil
, trapezoidRule
, trapezoidStencil
, simpsonsRule
, simpsonsStencil
, integrateWith
, integralApprox
, fromFunction
, fromFunctionMidpoint
) where
import Data.Coerce
import Data.Massiv.Array.Delayed.Pull (D)
import Data.Massiv.Array.Delayed.Windowed (DW)
import Data.Massiv.Array.Manifest.Internal
import Data.Massiv.Array.Ops.Construct (rangeInclusive)
import Data.Massiv.Array.Ops.Transform (extract')
import Data.Massiv.Array.Stencil
import Data.Massiv.Core.Common
midpointStencil ::
(Fractional e, Index ix)
=> e
-> Dim
-> Int
-> Stencil ix e e
midpointStencil :: e -> Dim -> Int -> Stencil ix e e
midpointStencil e
dx Dim
dim Int
k =
e -> Sz ix -> ix -> ((ix -> Value e) -> Value e) -> Stencil ix e e
forall ix e a.
Index ix =>
e -> Sz ix -> ix -> ((ix -> Value e) -> Value a) -> Stencil ix e a
makeStencilDef e
0 (ix -> Sz ix
forall ix. Index ix => ix -> Sz ix
Sz (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' (Int -> ix
forall ix. Index ix => Int -> ix
pureIndex Int
1) Dim
dim Int
k)) ix
forall ix. Index ix => ix
zeroIndex (((ix -> Value e) -> Value e) -> Stencil ix e e)
-> ((ix -> Value e) -> Value e) -> Stencil ix e e
forall a b. (a -> b) -> a -> b
$ \ix -> Value e
g ->
e -> Value e
forall (f :: * -> *) a. Applicative f => a -> f a
pure e
dx Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
* Int
-> (Int -> Bool)
-> (Int -> Int)
-> Value e
-> (Int -> Value e -> Value e)
-> Value e
forall a.
Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> a) -> a
loop Int
0 (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
k) (Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Value e
0 (\Int
i -> (Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
+ ix -> Value e
g (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim Int
i)))
{-# INLINE midpointStencil #-}
trapezoidStencil ::
(Fractional e, Index ix)
=> e
-> Dim
-> Int
-> Stencil ix e e
trapezoidStencil :: e -> Dim -> Int -> Stencil ix e e
trapezoidStencil e
dx Dim
dim Int
n =
e -> Sz ix -> ix -> ((ix -> Value e) -> Value e) -> Stencil ix e e
forall ix e a.
Index ix =>
e -> Sz ix -> ix -> ((ix -> Value e) -> Value a) -> Stencil ix e a
makeStencilDef e
0 (ix -> Sz ix
forall ix. Index ix => ix -> Sz ix
Sz (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' (Int -> ix
forall ix. Index ix => Int -> ix
pureIndex Int
1) Dim
dim (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))) ix
forall ix. Index ix => ix
zeroIndex (((ix -> Value e) -> Value e) -> Stencil ix e e)
-> ((ix -> Value e) -> Value e) -> Stencil ix e e
forall a b. (a -> b) -> a -> b
$ \ix -> Value e
g ->
e -> Value e
forall (f :: * -> *) a. Applicative f => a -> f a
pure e
dx Value e -> Value e -> Value e
forall a. Fractional a => a -> a -> a
/ Value e
2 Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
*
(Int
-> (Int -> Bool)
-> (Int -> Int)
-> Value e
-> (Int -> Value e -> Value e)
-> Value e
forall a.
Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> a) -> a
loop Int
1 (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n) (Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (ix -> Value e
g ix
forall ix. Index ix => ix
zeroIndex) (\Int
i -> (Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
+ Value e
2 Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
* ix -> Value e
g (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim Int
i))) Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
+
ix -> Value e
g (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim Int
n))
{-# INLINE trapezoidStencil #-}
simpsonsStencil ::
(Fractional e, Index ix)
=> e
-> Dim
-> Int
-> Stencil ix e e
simpsonsStencil :: e -> Dim -> Int -> Stencil ix e e
simpsonsStencil e
dx Dim
dim Int
n
| Int -> Bool
forall a. Integral a => a -> Bool
odd Int
n =
[Char] -> Stencil ix e e
forall a. HasCallStack => [Char] -> a
error ([Char] -> Stencil ix e e) -> [Char] -> Stencil ix e e
forall a b. (a -> b) -> a -> b
$
[Char]
"Number of sample points for Simpson's rule stencil should be even, but received: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
n
| Bool
otherwise =
e -> Sz ix -> ix -> ((ix -> Value e) -> Value e) -> Stencil ix e e
forall ix e a.
Index ix =>
e -> Sz ix -> ix -> ((ix -> Value e) -> Value a) -> Stencil ix e a
makeStencilDef e
0 (ix -> Sz ix
forall ix. Index ix => ix -> Sz ix
Sz (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' (Int -> ix
forall ix. Index ix => Int -> ix
pureIndex Int
1) Dim
dim (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1))) ix
forall ix. Index ix => ix
zeroIndex (((ix -> Value e) -> Value e) -> Stencil ix e e)
-> ((ix -> Value e) -> Value e) -> Stencil ix e e
forall a b. (a -> b) -> a -> b
$ \ix -> Value e
g ->
let simAcc :: Int -> (Value e, Value e) -> (Value e, Value e)
simAcc Int
i (Value e
prev, Value e
acc) =
let !fx3 :: Value e
fx3 = ix -> Value e
g (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2))
!newAcc :: Value e
newAcc = Value e
acc Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
+ Value e
prev Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
+ Value e
4 Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
* ix -> Value e
g (ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)) Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
+ Value e
fx3
in (Value e
fx3, Value e
newAcc)
in e -> Value e
forall (f :: * -> *) a. Applicative f => a -> f a
pure e
dx Value e -> Value e -> Value e
forall a. Fractional a => a -> a -> a
/ Value e
3 Value e -> Value e -> Value e
forall a. Num a => a -> a -> a
* (Value e, Value e) -> Value e
forall a b. (a, b) -> b
snd (Int
-> (Int -> Bool)
-> (Int -> Int)
-> (Value e, Value e)
-> (Int -> (Value e, Value e) -> (Value e, Value e))
-> (Value e, Value e)
forall a.
Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> a) -> a
loop Int
2 (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2) (Int -> (Value e, Value e) -> (Value e, Value e)
simAcc Int
0 (ix -> Value e
g ix
forall ix. Index ix => ix
zeroIndex, Value e
0)) Int -> (Value e, Value e) -> (Value e, Value e)
simAcc)
{-# INLINE simpsonsStencil #-}
integrateWith ::
(Fractional e, StrideLoad DW ix e, Mutable r ix e)
=> (Dim -> Int -> Stencil ix e e)
-> Dim
-> Int
-> Array r ix e
-> Array r ix e
integrateWith :: (Dim -> Int -> Stencil ix e e)
-> Dim -> Int -> Array r ix e -> Array r ix e
integrateWith Dim -> Int -> Stencil ix e e
stencil Dim
dim Int
n Array r ix e
arr =
Stride ix -> Array DW ix e -> Array r ix e
forall r ix e r'.
(Mutable r ix e, StrideLoad r' ix e) =>
Stride ix -> Array r' ix e -> Array r ix e
computeWithStride (ix -> Stride ix
forall ix. Index ix => ix -> Stride ix
Stride ix
nsz) (Array DW ix e -> Array r ix e) -> Array DW ix e -> Array r ix e
forall a b. (a -> b) -> a -> b
$ Border e -> Stencil ix e e -> Array r ix e -> Array DW ix e
forall r ix e a.
(Source r ix e, Manifest r ix e) =>
Border e -> Stencil ix e a -> Array r ix e -> Array DW ix a
mapStencil (e -> Border e
forall e. e -> Border e
Fill e
0) (Dim -> Int -> Stencil ix e e
stencil Dim
dim Int
n) Array r ix e
arr
where
!nsz :: ix
nsz = ix -> Dim -> Int -> ix
forall ix. Index ix => ix -> Dim -> Int -> ix
setDim' (Int -> ix
forall ix. Index ix => Int -> ix
pureIndex Int
1) Dim
dim Int
n
{-# INLINE integrateWith #-}
integralApprox ::
(Fractional e, StrideLoad DW ix e, Mutable r ix e)
=> (e -> Dim -> Int -> Stencil ix e e)
-> e
-> Sz ix
-> Int
-> Array r ix e
-> Array M ix e
integralApprox :: (e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
integralApprox e -> Dim -> Int -> Stencil ix e e
stencil e
d Sz ix
sz Int
n Array r ix e
arr =
ix -> Sz ix -> Array M ix e -> Array (R M) ix e
forall r ix e.
Extract r ix e =>
ix -> Sz ix -> Array r ix e -> Array (R r) ix e
extract' ix
forall ix. Index ix => ix
zeroIndex Sz ix
sz (Array M ix e -> Array (R M) ix e)
-> Array M ix e -> Array (R M) ix e
forall a b. (a -> b) -> a -> b
$ Array r ix e -> Array M ix e
forall r ix e. Manifest r ix e => Array r ix e -> Array M ix e
toManifest (Array r ix e -> Array M ix e) -> Array r ix e -> Array M ix e
forall a b. (a -> b) -> a -> b
$ Int
-> (Int -> Bool)
-> (Int -> Int)
-> Array r ix e
-> (Int -> Array r ix e -> Array r ix e)
-> Array r ix e
forall a.
Int -> (Int -> Bool) -> (Int -> Int) -> a -> (Int -> a -> a) -> a
loop Int
1 (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Dim -> Int
coerce (Sz ix -> Dim
forall ix (proxy :: * -> *). Index ix => proxy ix -> Dim
dimensions Sz ix
sz)) (Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Array r ix e
arr Int -> Array r ix e -> Array r ix e
forall r. Mutable r ix e => Int -> Array r ix e -> Array r ix e
integrateAlong
where
!dx :: e
dx = e
d e -> e -> e
forall a. Fractional a => a -> a -> a
/ Int -> e
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
integrateAlong :: Int -> Array r ix e -> Array r ix e
integrateAlong Int
dim = (Dim -> Int -> Stencil ix e e)
-> Dim -> Int -> Array r ix e -> Array r ix e
forall e ix r.
(Fractional e, StrideLoad DW ix e, Mutable r ix e) =>
(Dim -> Int -> Stencil ix e e)
-> Dim -> Int -> Array r ix e -> Array r ix e
integrateWith (e -> Dim -> Int -> Stencil ix e e
stencil e
dx) (Int -> Dim
Dim Int
dim) Int
n
{-# INLINE integrateAlong #-}
{-# INLINE integralApprox #-}
midpointRule ::
(Fractional e, StrideLoad DW ix e, Mutable r ix e)
=> Comp
-> r
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array M ix e
midpointRule :: Comp
-> r
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array M ix e
midpointRule Comp
comp r
r (Int -> e) -> ix -> e
f e
a e
d Sz ix
sz Int
n =
(e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
forall e ix r.
(Fractional e, StrideLoad DW ix e, Mutable r ix e) =>
(e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
integralApprox e -> Dim -> Int -> Stencil ix e e
forall e ix.
(Fractional e, Index ix) =>
e -> Dim -> Int -> Stencil ix e e
midpointStencil e
d Sz ix
sz Int
n (Array r ix e -> Array M ix e) -> Array r ix e -> Array M ix e
forall a b. (a -> b) -> a -> b
$ r -> Array D ix e -> Array r ix e
forall r ix e r'.
(Mutable r ix e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs r
r (Array D ix e -> Array r ix e) -> Array D ix e -> Array r ix e
forall a b. (a -> b) -> a -> b
$ Comp
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array D ix e
forall ix a e.
(Index ix, Fractional a) =>
Comp
-> ((Int -> a) -> ix -> e)
-> a
-> a
-> Sz ix
-> Int
-> Array D ix e
fromFunctionMidpoint Comp
comp (Int -> e) -> ix -> e
f e
a e
d Sz ix
sz Int
n
{-# INLINE midpointRule #-}
trapezoidRule ::
(Fractional e, StrideLoad DW ix e, Mutable r ix e)
=> Comp
-> r
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array M ix e
trapezoidRule :: Comp
-> r
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array M ix e
trapezoidRule Comp
comp r
r (Int -> e) -> ix -> e
f e
a e
d Sz ix
sz Int
n =
(e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
forall e ix r.
(Fractional e, StrideLoad DW ix e, Mutable r ix e) =>
(e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
integralApprox e -> Dim -> Int -> Stencil ix e e
forall e ix.
(Fractional e, Index ix) =>
e -> Dim -> Int -> Stencil ix e e
trapezoidStencil e
d Sz ix
sz Int
n (Array r ix e -> Array M ix e) -> Array r ix e -> Array M ix e
forall a b. (a -> b) -> a -> b
$ r -> Array D ix e -> Array r ix e
forall r ix e r'.
(Mutable r ix e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs r
r (Array D ix e -> Array r ix e) -> Array D ix e -> Array r ix e
forall a b. (a -> b) -> a -> b
$ Comp
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array D ix e
forall ix a e.
(Index ix, Fractional a) =>
Comp
-> ((Int -> a) -> ix -> e)
-> a
-> a
-> Sz ix
-> Int
-> Array D ix e
fromFunction Comp
comp (Int -> e) -> ix -> e
f e
a e
d Sz ix
sz Int
n
{-# INLINE trapezoidRule #-}
simpsonsRule ::
(Fractional e, StrideLoad DW ix e, Mutable r ix e)
=> Comp
-> r
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array M ix e
simpsonsRule :: Comp
-> r
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array M ix e
simpsonsRule Comp
comp r
r (Int -> e) -> ix -> e
f e
a e
d Sz ix
sz Int
n =
(e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
forall e ix r.
(Fractional e, StrideLoad DW ix e, Mutable r ix e) =>
(e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array M ix e
integralApprox e -> Dim -> Int -> Stencil ix e e
forall e ix.
(Fractional e, Index ix) =>
e -> Dim -> Int -> Stencil ix e e
simpsonsStencil e
d Sz ix
sz Int
n (Array r ix e -> Array M ix e) -> Array r ix e -> Array M ix e
forall a b. (a -> b) -> a -> b
$ r -> Array D ix e -> Array r ix e
forall r ix e r'.
(Mutable r ix e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs r
r (Array D ix e -> Array r ix e) -> Array D ix e -> Array r ix e
forall a b. (a -> b) -> a -> b
$ Comp
-> ((Int -> e) -> ix -> e)
-> e
-> e
-> Sz ix
-> Int
-> Array D ix e
forall ix a e.
(Index ix, Fractional a) =>
Comp
-> ((Int -> a) -> ix -> e)
-> a
-> a
-> Sz ix
-> Int
-> Array D ix e
fromFunction Comp
comp (Int -> e) -> ix -> e
f e
a e
d Sz ix
sz Int
n
{-# INLINE simpsonsRule #-}
fromFunction ::
(Index ix, Fractional a)
=> Comp
-> ((Int -> a) -> ix -> e)
-> a
-> a
-> Sz ix
-> Int
-> Array D ix e
fromFunction :: Comp
-> ((Int -> a) -> ix -> e)
-> a
-> a
-> Sz ix
-> Int
-> Array D ix e
fromFunction Comp
comp (Int -> a) -> ix -> e
f a
a a
d (Sz ix
sz) Int
n =
(Int -> a) -> ix -> e
f Int -> a
forall a. Integral a => a -> a
scale (ix -> e) -> Array D ix ix -> Array D ix e
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Comp -> ix -> ix -> Array D ix ix
forall ix. Index ix => Comp -> ix -> ix -> Array D ix ix
rangeInclusive Comp
comp ix
forall ix. Index ix => ix
zeroIndex ((Int -> Int) -> ix -> ix
forall ix. Index ix => (Int -> Int) -> ix -> ix
liftIndex (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
*) ix
sz)
where
nFrac :: a
nFrac = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
scale :: a -> a
scale a
i = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
d a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
nFrac
{-# INLINE scale #-}
{-# INLINE fromFunction #-}
fromFunctionMidpoint
:: (Index ix, Fractional a) =>
Comp -> ((Int -> a) -> ix -> e) -> a -> a -> Sz ix -> Int -> Array D ix e
fromFunctionMidpoint :: Comp
-> ((Int -> a) -> ix -> e)
-> a
-> a
-> Sz ix
-> Int
-> Array D ix e
fromFunctionMidpoint Comp
comp (Int -> a) -> ix -> e
f a
a a
d (Sz ix
sz) Int
n =
(Int -> a) -> ix -> e
f Int -> a
forall a. Integral a => a -> a
scale (ix -> e) -> Array D ix ix -> Array D ix e
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Comp -> ix -> ix -> Array D ix ix
forall ix. Index ix => Comp -> ix -> ix -> Array D ix ix
rangeInclusive Comp
comp ix
forall ix. Index ix => ix
zeroIndex ((Int -> Int) -> ix -> ix
forall ix. Index ix => (Int -> Int) -> ix -> ix
liftIndex (\Int
i -> Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) ix
sz)
where
nFrac :: a
nFrac = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
dx2 :: a
dx2 = a
d a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
nFrac a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2
scale :: a -> a
scale a
i = a
dx2 a -> a -> a
forall a. Num a => a -> a -> a
+ a
a a -> a -> a
forall a. Num a => a -> a -> a
+ a
d a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
nFrac
{-# INLINE scale #-}
{-# INLINE fromFunctionMidpoint #-}