{-# LANGUAGE BangPatterns     #-}
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module      : Data.Massiv.Array.Numeric.Integral
-- Copyright   : (c) Alexey Kuleshevich 2018-2019
-- License     : BSD3
-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
-- Stability   : experimental
-- Portability : non-portable
--
module Data.Massiv.Array.Numeric.Integral
  (
  -- $integral_intro
  --
  -- * Integral Approximation
  -- ** Midpoint Rule
    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
  -- * Helper functions
  ) 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


-- |
--
-- __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 dx dim k =
  makeStencilDef 0 (Sz (setDim' (pureIndex 1) dim k)) zeroIndex $ \g ->
    pure dx * loop 0 (< k) (+ 1) 0 (\i -> (+ g (setDim' zeroIndex dim 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 dx dim n =
  makeStencilDef 0 (Sz (setDim' (pureIndex 1) dim (n + 1))) zeroIndex $ \g ->
    pure dx / 2 *
    (loop 1 (< n) (+ 1) (g zeroIndex) (\i -> (+ 2 * g (setDim' zeroIndex dim i))) +
     g (setDim' zeroIndex dim 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 dx dim n
  | odd n =
    error $
    "Number of sample points for Simpson's rule stencil should be even, but received: " ++ show n
  | otherwise =
    makeStencilDef 0 (Sz (setDim' (pureIndex 1) dim (n + 1))) zeroIndex $ \g ->
      let simAcc i (prev, acc) =
            let !fx3 = g (setDim' zeroIndex dim (i + 2))
                !newAcc = acc + prev + 4 * g (setDim' zeroIndex dim (i + 1)) + fx3
             in (fx3, newAcc)
       in pure dx / 3 * snd (loop 2 (< n - 1) (+ 2) (simAcc 0 (g zeroIndex, 0)) simAcc)
{-# INLINE simpsonsStencil #-}

-- | Integrate with a stencil along a particular dimension.
integrateWith ::
     (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
integrateWith stencil dim n arr =
  computeWithStride (Stride nsz) $ mapStencil (Fill 0) (stencil dim n) arr
  where
    !nsz = setDim' (pureIndex 1) dim n
{-# INLINE integrateWith #-}


-- | Compute an approximation of integral using a supplied rule in a form of `Stencil`.
integralApprox ::
     (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
integralApprox stencil d sz n arr =
  extract' zeroIndex sz $ toManifest $ loop 1 (<= coerce (dimensions sz)) (+ 1) arr applyStencil
  where
    !dx = d / fromIntegral n
    applyStencil dim = integrateWith (stencil dx) (Dim dim) n
    {-# INLINE applyStencil #-}
{-# INLINE integralApprox #-}


-- | Use midpoint rule to approximate an integral.
midpointRule ::
     (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
midpointRule comp r f a d sz n =
  integralApprox midpointStencil d sz n $ computeAs r $ fromFunctionMidpoint comp f a d sz n
{-# INLINE midpointRule #-}


-- | Use trapezoid rule to approximate an integral.
trapezoidRule ::
     (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
trapezoidRule comp r f a d sz n =
  integralApprox trapezoidStencil d sz n $ computeAs r $ fromFunction comp f a d sz n
{-# INLINE trapezoidRule #-}

-- | Use Simpson's rule to approximate an integral.
simpsonsRule ::
     (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
simpsonsRule comp r f a d sz n =
  integralApprox simpsonsStencil d sz n $ computeAs r $ fromFunction comp f a d sz 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 comp f a d (Sz sz) n =
  f scale <$> rangeInclusive comp zeroIndex (liftIndex (n *) sz)
  where
    nFrac = fromIntegral n
    scale i = a + d * fromIntegral i / 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 comp f a d (Sz sz) n =
  f scale <$> rangeInclusive comp zeroIndex (liftIndex (\i -> n * i - 1) sz)
  where
    nFrac = fromIntegral n
    dx2 = d / nFrac / 2
    scale i = dx2 + a + d * fromIntegral i / 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, 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 ]