{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}

-- |
-- Module      : Data.Massiv.Array.Numeric.Integral
-- Copyright   : (c) Alexey Kuleshevich 2018-2022
-- License     : BSD3
-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
-- Stability   : experimental
-- Portability : non-portable
module Data.Massiv.Array.Numeric.Integral (
  -- $integral_intro
  midpointRule,
  midpointStencil,

  -- ** Trapezoid Rule
  trapezoidRule,
  trapezoidStencil,

  -- ** Simpson's Rule
  simpsonsRule,
  simpsonsStencil,

  -- * General Integral approximation
  integrateWith,
  integralApprox,

  -- * From functions

  -- ** Sampled at the edge
  fromFunction,

  -- ** Sampled at the midpoint
  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

-- |
--
-- __Midpoint Rule__
--
-- \[
-- \int_{{\,a}}^{{\,b}}{{f\left( x \right)\,dx}} \approx \Delta x \cdot \,f\left( {x_1 + \frac{\Delta x}{2}} \right) + \Delta x \cdot \,f\left( {x_2 + \frac{\Delta x}{2}} \right) +  \cdots  + \Delta x \cdot \,f\left( {x_n + \frac{\Delta x}{2}} \right)
-- \]
midpointStencil
  :: (Fractional e, Index ix)
  => e
  -- ^ @Δx@ - distance between sample points
  -> Dim
  -- ^ Dimension along which to integrate
  -> Int
  -- ^ @n@ - number of sample points.
  -> 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 =
  Sz ix -> ix -> (ix -> (ix -> e) -> e) -> Stencil ix e e
forall ix e a.
Index ix =>
Sz ix -> ix -> (ix -> (ix -> e) -> a) -> Stencil ix e a
makeUnsafeStencil (ix -> Sz ix
forall ix. Index ix => ix -> Sz ix
Sz (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, 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 -> (ix -> e) -> e) -> Stencil ix e e)
-> (ix -> (ix -> e) -> e) -> Stencil ix e e
forall a b. (a -> b) -> a -> b
$ \ix
_ ix -> e
g ->
    e
dx e -> e -> e
forall a. Num a => a -> a -> a
* Int -> (Int -> Bool) -> (Int -> Int) -> e -> (Int -> e -> e) -> 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) e
0 (\Int
i -> (e -> e -> e
forall a. Num a => a -> a -> a
+ ix -> e
g (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, Index ix) => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim Int
i)))
{-# INLINE midpointStencil #-}

-- |
--
-- __Trapezoid Rule__
--
-- \[
-- \int_{{\,a}}^{{\,b}}{{f\left( x \right)\,dx}} \approx \frac{{\Delta x}}{2}\cdot\left( {f\left( {{x_0}} \right) + f\left( {{x_1}} \right)} \right) + \frac{{\Delta x}}{2}\cdot\left( {f\left( {{x_1}} \right) + f\left( {{x_2}} \right)} \right) +  \cdots  + \frac{{\Delta x}}{2}\cdot\left( {f\left( {{x_{n - 1}}} \right) + f\left( {{x_n}} \right)} \right)
-- \]
trapezoidStencil
  :: (Fractional e, Index ix)
  => e
  -- ^ @Δx@ - distance between sample points
  -> Dim
  -- ^ Dimension along which to integrate
  -> Int
  -- ^ @n@ - number of sample points.
  -> 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 =
  Sz ix -> ix -> (ix -> (ix -> e) -> e) -> Stencil ix e e
forall ix e a.
Index ix =>
Sz ix -> ix -> (ix -> (ix -> e) -> a) -> Stencil ix e a
makeUnsafeStencil (ix -> Sz ix
forall ix. Index ix => ix -> Sz ix
Sz (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, 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 -> (ix -> e) -> e) -> Stencil ix e e)
-> (ix -> (ix -> e) -> e) -> Stencil ix e e
forall a b. (a -> b) -> a -> b
$ \ix
_ ix -> e
g ->
    (e
dx e -> e -> e
forall a. Fractional a => a -> a -> a
/ e
2)
      e -> e -> e
forall a. Num a => a -> a -> a
* ( Int -> (Int -> Bool) -> (Int -> Int) -> e -> (Int -> e -> e) -> 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 -> e
g ix
forall ix. Index ix => ix
zeroIndex) (\Int
i -> (e -> e -> e
forall a. Num a => a -> a -> a
+ e
2 e -> e -> e
forall a. Num a => a -> a -> a
* ix -> e
g (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, Index ix) => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim Int
i)))
            e -> e -> e
forall a. Num a => a -> a -> a
+ ix -> e
g (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, Index ix) => ix -> Dim -> Int -> ix
setDim' ix
forall ix. Index ix => ix
zeroIndex Dim
dim Int
n)
        )
{-# INLINE trapezoidStencil #-}

-- |
--
-- __Simpson's Rule__
--
-- \[
-- \int_{{\,a}}^{{\,b}}{{f\left( x \right)\,dx}} \approx \frac{{\Delta x}}{3}\cdot\left( {f\left( {{x_0}} \right) + 4\cdot f\left( {{x_1}} \right) + f\left( {{x_2}} \right)} \right) + \frac{{\Delta x}}{3}\cdot\left( {f\left( {{x_2}} \right) + 4\cdot f\left( {{x_3}} \right) + f\left( {{x_4}} \right)} \right) +  \cdots  + \frac{{\Delta x}}{3}\cdot\left( {f\left( {{x_{n - 2}}} \right) + 4\cdot f\left( {{x_{n - 1}}} \right) + f\left( {{x_n}} \right)} \right)
-- \]
simpsonsStencil
  :: (Fractional e, Index ix)
  => e
  -- ^ @Δx@ - distance between sample points
  -> Dim
  -- ^ Dimension along which to integrate
  -> Int
  -- ^ @n@ - Number of sample points. This value should be even, otherwise error.
  -> 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
  | 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 =
      Sz ix -> ix -> (ix -> (ix -> e) -> e) -> Stencil ix e e
forall ix e a.
Index ix =>
Sz ix -> ix -> (ix -> (ix -> e) -> a) -> Stencil ix e a
makeUnsafeStencil (ix -> Sz ix
forall ix. Index ix => ix -> Sz ix
Sz (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, 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 -> (ix -> e) -> e) -> Stencil ix e e)
-> (ix -> (ix -> e) -> e) -> Stencil ix e e
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 (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, 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 :: e
newAcc = e
acc e -> e -> e
forall a. Num a => a -> a -> a
+ e
prev e -> e -> e
forall a. Num a => a -> a -> a
+ e
4 e -> e -> e
forall a. Num a => a -> a -> a
* ix -> e
g (ix -> Dim -> Int -> ix
forall ix. (HasCallStack, 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)) e -> e -> e
forall a. Num a => a -> a -> a
+ e
fx3
               in (e
fx3, e
newAcc)
         in e
dx e -> e -> e
forall a. Fractional a => a -> a -> a
/ e
3 e -> e -> e
forall a. Num a => a -> a -> a
* (e, e) -> e
forall a b. (a, b) -> b
snd (Int
-> (Int -> Bool)
-> (Int -> Int)
-> (e, e)
-> (Int -> (e, e) -> (e, e))
-> (e, 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 -> (e, e) -> (e, e)
simAcc Int
0 (ix -> e
g ix
forall ix. Index ix => ix
zeroIndex, e
0)) Int -> (e, e) -> (e, e)
simAcc)
{-# INLINE simpsonsStencil #-}

-- | Integrate with a stencil along a particular dimension.
integrateWith
  :: (Fractional e, StrideLoad DW ix e, Manifest r e)
  => (Dim -> Int -> Stencil ix e e)
  -> Dim
  -- ^ Dimension along which integration should be estimated.
  -> Int
  -- ^ @n@ - Number of samples
  -> 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 =
  Stride ix -> Array DW ix e -> Array r ix e
forall r ix e r'.
(Manifest r 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 ix r e a.
(Index ix, Manifest r 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. (HasCallStack, Index ix) => ix -> Dim -> Int -> ix
setDim' (Int -> ix
forall ix. Index ix => Int -> ix
pureIndex Int
1) Dim
dim Int
n
{-# INLINE integrateWith #-}

-- | Compute an approximation of integral using a supplied rule in a form of `Stencil`.
integralApprox
  :: (Fractional e, StrideLoad DW ix e, Manifest r e)
  => (e -> Dim -> Int -> Stencil ix e e)
  -- ^ Integration Stencil
  -> e
  -- ^ @d@ - Length of interval per cell
  -> Sz ix
  -- ^ @sz@ - Result size of the matrix
  -> Int
  -- ^ @n@ - Number of samples
  -> Array r ix e
  -- ^ Array with values of @f(x,y,..)@ that will be used as source for integration.
  -> 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 =
  ix -> Sz ix -> Array r ix e -> Array D ix e
forall r ix e.
(HasCallStack, Index ix, Source r e) =>
ix -> Sz ix -> Array r ix e -> Array D ix e
extract' ix
forall ix. Index ix => ix
zeroIndex Sz ix
sz (Array r ix e -> Array D ix e) -> Array r ix e -> Array D 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
forall a b. Coercible a b => a -> b
coerce (Sz ix -> Dim
forall ix (proxy :: * -> *). Index ix => proxy ix -> Dim
forall (proxy :: * -> *). 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}. Manifest r 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, 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 #-}

-- | Use midpoint rule to approximate an integral.
midpointRule
  :: (Fractional e, StrideLoad DW ix e, Manifest r e)
  => Comp
  -- ^ Computation strategy.
  -> r
  -- ^ Intermediate array representation.
  -> ((Int -> e) -> ix -> e)
  -- ^ @f(x,y,...)@ - Function to integrate
  -> e
  -- ^ @a@ - Starting value point.
  -> e
  -- ^ @d@ - Distance per matrix cell.
  -> Sz ix
  -- ^ @sz@ - Result matrix size.
  -> Int
  -- ^ @n@ - Number of sample points per cell in each direction.
  -> 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 =
  (e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array D ix e
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
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 D ix e) -> Array r ix e -> Array D ix e
forall a b. (a -> b) -> a -> b
$ r -> Array D ix e -> Array r ix e
forall r e r' ix.
(Manifest r 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 #-}

-- | Use trapezoid rule to approximate an integral.
trapezoidRule
  :: (Fractional e, StrideLoad DW ix e, Manifest r e)
  => Comp
  -- ^ Computation strategy
  -> r
  -- ^ Intermediate array representation
  -> ((Int -> e) -> ix -> e)
  -- ^ @f(x,y,...)@ - function to integrate
  -> e
  -- ^ @a@ - Starting value point.
  -> e
  -- ^ @d@ - Distance per matrix cell.
  -> Sz ix
  -- ^ @sz@ - Result matrix size.
  -> Int
  -- ^ @n@ - Number of sample points per cell in each direction.
  -> 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 =
  (e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array D ix e
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
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 D ix e) -> Array r ix e -> Array D ix e
forall a b. (a -> b) -> a -> b
$ r -> Array D ix e -> Array r ix e
forall r e r' ix.
(Manifest r 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 #-}

-- | Use Simpson's rule to approximate an integral.
simpsonsRule
  :: (Fractional e, StrideLoad DW ix e, Manifest r e)
  => Comp
  -- ^ Computation strategy
  -> r
  -- ^ Intermediate array representation
  -> ((Int -> e) -> ix -> e)
  -- ^ @f(x,y,...)@ - Function to integrate
  -> e
  -- ^ @a@ - Starting value point.
  -> e
  -- ^ @d@ - Distance per matrix cell.
  -> Sz ix
  -- ^ @sz@ - Result matrix size.
  -> Int
  -- ^ @n@ - Number of sample points per cell in each direction. This value must be even,
  -- otherwise error.
  -> 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 =
  (e -> Dim -> Int -> Stencil ix e e)
-> e -> Sz ix -> Int -> Array r ix e -> Array D ix e
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
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 D ix e) -> Array r ix e -> Array D ix e
forall a b. (a -> b) -> a -> b
$ r -> Array D ix e -> Array r ix e
forall r e r' ix.
(Manifest r 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 #-}

-- | Create an array from a function with sample points at the edges
--
-- >>> fromFunction Seq (\ scale (i :. j) -> scale i + scale j :: Double) (-2) 1 (Sz 4) 2
-- Array D Seq (Sz (9 :. 9))
--   [ [ -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0 ]
--   , [ -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5 ]
--   , [ -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0 ]
--   , [ -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5 ]
--   , [ -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0 ]
--   , [ -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 ]
--   , [ -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 ]
--   , [ -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 ]
--   , [ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 ]
--   ]
fromFunction
  :: (Index ix, Fractional a)
  => Comp
  -- ^ Computation strategy
  -> ((Int -> a) -> ix -> e)
  -- ^ A function that will produce elements of scaled up array. First argument is a scaling
  -- function that should be applied to individual indicies.
  -> a
  -- ^ @a@ - Starting point
  -> a
  -- ^ @d@ - Distance per cell
  -> Sz ix
  -- ^ @sz@ - Size of the desired array
  -> Int
  -- ^ @n@ - Scaling factor, i.e. number of sample points per cell.
  -> 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 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 #-}

-- | Similar to `fromFunction`, but will create an array from a function with sample points in the
-- middle of cells.
--
-- >>> fromFunctionMidpoint Seq (\ scale (i :. j) -> scale i + scale j :: Double) (-2) 1 (Sz 4) 2
-- Array D Seq (Sz (8 :. 8))
--   [ [ -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0 ]
--   , [ -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5 ]
--   , [ -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0 ]
--   , [ -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5 ]
--   , [ -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0 ]
--   , [ -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 ]
--   , [ -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 ]
--   , [ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 ]
--   ]
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 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 #-}

-- $integral_intro
--
-- Inspiration for the code in this module was taken from [Paul Dawkins Online
-- Notes](http://tutorial.math.lamar.edu). In particular the page on [Integral
-- Approximation](http://tutorial.math.lamar.edu/Classes/CalcII/ApproximatingDefIntegrals.aspx),
-- so if you need to brush up on some theory it is a great place to start.
--
-- Implementation-wise, integral approximation here relies heavily on stencils
-- with stride, because such computation is fast and is automatically
-- parallelizable.
--
-- Here are some examples of where this can be useful:
--
-- === Integral of a function on a region
--
-- Say we have a gaussian @f(x) = e^(x^2)@ on interval @[0, 2]@ (as in Paul's tutorial above). For
-- this we define a function @f@, an array with equally spaced (/dx/) sample input values and apply
-- the function to that array, which will give us an array of @n + 1@ sample points, or looking from
-- a different angle @n@ intervals.
--
-- >>> import Data.Massiv.Array
-- >>> f x = exp ( x ^ (2 :: Int) ) :: Float
-- >>> fromFunction Seq (\ scale x -> f (scale x)) 0 2 (Sz1 1) 4
-- Array D Seq (Sz1 5)
--   [ 1.0, 1.2840254, 2.7182817, 9.487736, 54.59815 ]
--
-- Once we have that array of sample points ready, we could use `integralApprox` and one of the
-- stencils to compute an integral, but there are already functions that will do both steps for you:
--
-- >>> simpsonsRule Seq U (\ scale x -> f (scale x)) 0 2 (Sz1 1) 4
-- Array D Seq (Sz1 1)
--   [ 17.353626 ]
--
-- @scale@ is the function that will change an array index into equally spaced and
-- appropriately shifted values of @x, y, ...@ before they can get applied to @f(x, y, ...)@
--
-- === Accurate values of a function
--
-- Another very useful place where integral approximation can be used is when a more accurate
-- representation of a non-linear function is desired. Consider the same gaussian function applied
-- to equally spaced values, with zero being in the middle of the vector:
--
-- >>> xArr = makeArrayR D Seq (Sz1 4) $ \ i -> fromIntegral i - 1.5 :: Float
-- >>> xArr
-- Array D Seq (Sz1 4)
--   [ -1.5, -0.5, 0.5, 1.5 ]
-- >>> fmap f xArr
-- Array D Seq (Sz1 4)
--   [ 9.487736, 1.2840254, 1.2840254, 9.487736 ]
--
-- The problem with above example is that computed values do not accurately represent the total
-- value contained within each vector cell. For that reason if your were to later use it for example
-- as convolution stencil, approximation would be very poor. The way to solve it is to approximate
-- an integral across each cell of vector by drastically blowing up the @xArr@ and then reducing it
-- to a smaller array by using one of the approximation rules:
--
-- >>> startValue = -2 :: Float
-- >>> distPerCell = 1 :: Float
-- >>> desiredSize = Sz1 4 :: Sz1
-- >>> numSamples = 4 :: Int
-- >>> xArrX4 = fromFunction Seq ($) startValue distPerCell desiredSize numSamples
-- >>> xArrX4
-- Array D Seq (Sz1 17)
--   [ -2.0, -1.75, -1.5, -1.25, -1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0 ]
-- >>> yArrX4 = computeAs U $ fmap f xArrX4
-- >>> integralApprox trapezoidStencil distPerCell desiredSize numSamples yArrX4
-- Array D Seq (Sz1 4)
--   [ 16.074406, 1.4906789, 1.4906789, 16.074408 ]
--
-- We can clearly see the difference is huge, but it doesn't mean it is much better than our
-- previous estimate. In order to get more accurate results we can use a better Simpson's rule for
-- approximation and many more sample points. There is no need to create individual arrays @xArrX4@
-- and @yArrX4@, there are functions like `simpsonsRule` that will take care of it for us:
--
-- >>> simpsonsRule Seq U (\ scale i -> f (scale i)) startValue distPerCell desiredSize 128
-- Array D Seq (Sz1 4)
--   [ 14.989977, 1.4626511, 1.4626517, 14.989977 ]