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 example 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.

>>> f x = exp ( x ^ 2 )
>>> fromFunction Seq (\ scale x -> f (scale x)) 0 2 1 4
(Array D Seq (5)
  [ 1.0,1.2840254166877414,2.718281828459045,9.487735836358526,54.598150033144236 ])

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 1 4
  [ 17.353626450374566 ])

Scaling function scale is what 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 (Ix1 4) $ \ i -> (fromIntegral i - 1.5 :: Float)
>>> xArr
(Array D Seq (4)
  [ -1.5,-0.5,0.5,1.5 ])
>>> fmap f xArr
(Array D Seq (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 = 4 :: Ix1
>>> numSamples = 4 :: Int
>>> xArrX4 = fromFunction Seq ($) startValue distPerCell desiredSize numSamples
>>> xArrX4
(Array D Seq (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 (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 (4)
  [ 14.989977,1.4626511,1.4626517,14.989977 ])

midpointRule

Arguments

:: (Fractional e, Load 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.

-> 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, Load 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 point

-> e

d - distance per matrix cell

-> ix

sz - end 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, Load 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 point

-> e

d - distance per matrix cell

-> ix

sz - end 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, Load 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, Load DW ix e, Mutable r ix e)

=> (e -> Dim -> Int -> Stencil ix e e)

Integration Stencil

-> e

d - Length of interval per cell

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

-> 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) (-2) 1 (4 :. 4) 2
(Array D Seq (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 -> 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) (-2) 1 (4 :. 4) 2
(Array D Seq (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 ]
  ])