Data.Massiv.Array.Numeric.Integral

Documentation

Inspiration for the code in this module was taken from Paul Dawkins Online Notes. In particular the page on Integral Approximation, 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, as 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 M 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 M 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 xArr and yArr, there are functions like simpsonRule that will take care it for you:

>>> simpsonsRule Seq U (\ scale i -> f (scale i)) startValue distPerCell desiredSize 128
Array M Seq (Sz1 4)
  [ 14.989977, 1.4626511, 1.4626517, 14.989977 ]

midpointRule

Arguments

:: (Fractional e, StrideLoad DW ix e, Mutable r ix 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 M ix e

Use midpoint rule to approximate an integral.

midpointStencil

Arguments

:: (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

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) \]

Trapezoid Rule

trapezoidRule

Arguments

:: (Fractional e, StrideLoad DW ix e, Mutable r ix 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 M ix e

Use trapezoid rule to approximate an integral.

trapezoidStencil

Arguments

:: (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

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) \]

Simpson's Rule

simpsonsRule

Arguments

:: (Fractional e, StrideLoad DW ix e, Mutable r ix 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 M ix e

Use Simpson's rule to approximate an integral.

simpsonsStencil

Arguments

:: (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

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) \]

General Integral approximation

integrateWith

Arguments

:: (Fractional e, StrideLoad DW ix e, Mutable r ix 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

Integrate with a stencil along a particular dimension.

integralApprox

Arguments

:: (Fractional e, StrideLoad DW ix e, Mutable r ix 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 M ix e

Compute an approximation of integral using a supplied rule in a form of Stencil.

From functions

Sampled at the edge

fromFunction

Arguments

:: (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

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 ]
  ]

Sampled at the midpoint

fromFunctionMidpoint :: (Index ix, Fractional a) => Comp -> ((Int -> a) -> ix -> e) -> a -> a -> Sz ix -> Int -> Array D ix e

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 ]
  ]

Helper functions