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