{-# LANGUAGE GeneralizedNewtypeDeriving #-}

{-|
Copyright   : Predictable Network Solutions Ltd., 2020-2024
License     : BSD-3-Clause
Description : Finite signed measures on the number line.
-}
module Numeric.Measure.Finite.Mixed
    ( -- * Type
      Measure
    , zero
    , dirac
    , uniform
    , distribution
    , fromDistribution

      -- * Observations
    , total
    , support
    , isPositive
    , integrate

      -- * Operations, numerical
    , add
    , scale
    , translate
    , convolve
    ) where

import Data.Function.Class
    ( Function (..)
    )
import Data.List
    ( scanl'
    )
import Control.DeepSeq
    ( NFData
    )
import Numeric.Function.Piecewise
    ( Piecewise
    )
import Numeric.Polynomial.Simple
    ( Poly
    )

import qualified Data.Map.Strict as Map
import qualified Numeric.Function.Piecewise as Piecewise
import qualified Numeric.Measure.Discrete as D
import qualified Numeric.Polynomial.Simple as Poly

{-----------------------------------------------------------------------------
    Type
------------------------------------------------------------------------------}
-- | A finite
-- [signed measure](https://en.wikipedia.org/wiki/Signed_measure)
-- on the number line.
newtype Measure a = Measure (Piecewise (Poly a))
    -- INVARIANT: Adjacent pieces contain distinct objects.
    -- INVARIANT: The last piece is a constant polynomial,
    --            so that the measure is finite.
    deriving (Int -> Measure a -> ShowS
[Measure a] -> ShowS
Measure a -> String
(Int -> Measure a -> ShowS)
-> (Measure a -> String)
-> ([Measure a] -> ShowS)
-> Show (Measure a)
forall a. Show a => Int -> Measure a -> ShowS
forall a. Show a => [Measure a] -> ShowS
forall a. Show a => Measure a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall a. Show a => Int -> Measure a -> ShowS
showsPrec :: Int -> Measure a -> ShowS
$cshow :: forall a. Show a => Measure a -> String
show :: Measure a -> String
$cshowList :: forall a. Show a => [Measure a] -> ShowS
showList :: [Measure a] -> ShowS
Show, Measure a -> ()
(Measure a -> ()) -> NFData (Measure a)
forall a. NFData a => Measure a -> ()
forall a. (a -> ()) -> NFData a
$crnf :: forall a. NFData a => Measure a -> ()
rnf :: Measure a -> ()
NFData)

-- | @eval (distribution m) x@ is the measure of the interval \( (-∞, x] \).
--
-- This is known as the [distribution function
-- ](https://en.wikipedia.org/wiki/Distribution_function_%28measure_theory%29).
distribution :: (Ord a, Num a) => Measure a -> Piecewise (Poly a)
distribution :: forall a. (Ord a, Num a) => Measure a -> Piecewise (Poly a)
distribution (Measure Piecewise (Poly a)
p) = Piecewise (Poly a)
p

-- | Construct a signed measure from its
-- [distribution function
-- ](https://en.wikipedia.org/wiki/Distribution_function_%28measure_theory%29).
--
-- Return 'Nothing' if the measure is not finite,
-- that is if the last piece of the piecewise function is not constant.
fromDistribution
    :: (Ord a, Num a)
    => Piecewise (Poly a) -> Maybe (Measure a)
fromDistribution :: forall a. (Ord a, Num a) => Piecewise (Poly a) -> Maybe (Measure a)
fromDistribution Piecewise (Poly a)
pieces
    | Piecewise (Poly a) -> Bool
forall a. (Ord a, Num a) => Piecewise (Poly a) -> Bool
isEventuallyConstant Piecewise (Poly a)
pieces = Measure a -> Maybe (Measure a)
forall a. a -> Maybe a
Just (Measure a -> Maybe (Measure a)) -> Measure a -> Maybe (Measure a)
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a) -> Piecewise (Poly a)
forall a.
(Ord a, Num a) =>
Piecewise (Poly a) -> Piecewise (Poly a)
trim Piecewise (Poly a)
pieces
    | Bool
otherwise = Maybe (Measure a)
forall a. Maybe a
Nothing

-- | Test whether a piecewise polynomial is consant as x -> ∞.
isEventuallyConstant :: (Ord a, Num a) => Piecewise (Poly a) -> Bool
isEventuallyConstant :: forall a. (Ord a, Num a) => Piecewise (Poly a) -> Bool
isEventuallyConstant Piecewise (Poly a)
pieces
    | [(a, Poly a)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, Poly a)]
xpolys = Bool
True
    | Bool
otherwise = (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0) (Int -> Bool) -> ((a, Poly a) -> Int) -> (a, Poly a) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Poly a -> Int
forall a. (Eq a, Num a) => Poly a -> Int
Poly.degree (Poly a -> Int) -> ((a, Poly a) -> Poly a) -> (a, Poly a) -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a, Poly a) -> Poly a
forall a b. (a, b) -> b
snd ((a, Poly a) -> Bool) -> (a, Poly a) -> Bool
forall a b. (a -> b) -> a -> b
$ [(a, Poly a)] -> (a, Poly a)
forall a. HasCallStack => [a] -> a
last [(a, Poly a)]
xpolys
  where
    xpolys :: [(Domain (Poly a), Poly a)]
xpolys = Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
pieces

-- | Internal.
-- Join all intervals whose polynomials are equal.
trim :: (Ord a, Num a) => Piecewise (Poly a) -> Piecewise (Poly a)
trim :: forall a.
(Ord a, Num a) =>
Piecewise (Poly a) -> Piecewise (Poly a)
trim = Piecewise (Poly a) -> Piecewise (Poly a)
forall o. (Eq o, Num o) => Piecewise o -> Piecewise o
Piecewise.trim

-- | Two measures are equal if they yield the same measures on every set.
--
-- > mx == my
-- >   implies
-- >   forall t. eval (distribution mx) t = eval (distribution my) t
instance (Ord a, Num a) => Eq (Measure a) where
    Measure Piecewise (Poly a)
mx == :: Measure a -> Measure a -> Bool
== Measure Piecewise (Poly a)
my =
        Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
mx [(a, Poly a)] -> [(a, Poly a)] -> Bool
forall a. Eq a => a -> a -> Bool
== Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
my

{-----------------------------------------------------------------------------
    Operations
------------------------------------------------------------------------------}
-- | The measure that assigns @0@ to every set.
zero :: Num a => Measure a
zero :: forall a. Num a => Measure a
zero = Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure Piecewise (Poly a)
forall o. Piecewise o
Piecewise.zero

-- | A
-- [Dirac measure](https://en.wikipedia.org/wiki/Dirac_measure)
-- at the given point @x@.
--
-- > total (dirac x) = 1
dirac :: (Ord a, Num a) => a -> Measure a
dirac :: forall a. (Ord a, Num a) => a -> Measure a
dirac a
x = Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ [(Domain (Poly a), Poly a)] -> Piecewise (Poly a)
forall o. Ord (Domain o) => [(Domain o, o)] -> Piecewise o
Piecewise.fromAscPieces [(a
Domain (Poly a)
x, a -> Poly a
forall a. a -> Poly a
Poly.constant a
1)]

-- | The probability measure of a uniform probability distribution
-- in the interval \( [x,y) \).
--
-- > total (uniform x y) = 1
uniform :: (Ord a, Num a, Fractional a) => a -> a -> Measure a
uniform :: forall a. (Ord a, Num a, Fractional a) => a -> a -> Measure a
uniform a
x a
y = Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
x a
y of
    Ordering
EQ -> [(Domain (Poly a), Poly a)] -> Piecewise (Poly a)
forall o. Ord (Domain o) => [(Domain o, o)] -> Piecewise o
Piecewise.fromAscPieces [(a
Domain (Poly a)
x, Poly a
1)]
    Ordering
_  -> [(Domain (Poly a), Poly a)] -> Piecewise (Poly a)
forall o. Ord (Domain o) => [(Domain o, o)] -> Piecewise o
Piecewise.fromAscPieces [(a
Domain (Poly a)
low, Poly a
poly), (a
Domain (Poly a)
high, Poly a
1)]
  where
    low :: a
low = a -> a -> a
forall a. Ord a => a -> a -> a
min a
x a
y
    high :: a
high = a -> a -> a
forall a. Ord a => a -> a -> a
max a
x a
y
    poly :: Poly a
poly = (a, a) -> (a, a) -> Poly a
forall a. (Eq a, Fractional a) => (a, a) -> (a, a) -> Poly a
Poly.lineFromTo (a
low, a
0) (a
high, a
1)

-- | The total of the measure applied to the set of real numbers.
total :: (Ord a, Num a) => Measure a -> a
total :: forall a. (Ord a, Num a) => Measure a -> a
total (Measure Piecewise (Poly a)
p) =
    case Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
p of
        [] -> a
0
        [(Domain (Poly a), Poly a)]
ps -> Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval ((a, Poly a) -> Poly a
forall a b. (a, b) -> b
snd ([(a, Poly a)] -> (a, Poly a)
forall a. HasCallStack => [a] -> a
last [(a, Poly a)]
[(Domain (Poly a), Poly a)]
ps)) a
Domain (Poly a)
0

-- | The 'support' is the smallest closed, contiguous interval \( [x,y] \)
-- outside of which the measure is zero.
--
-- Returns 'Nothing' if the interval is empty.
support :: (Ord a, Num a) => Measure a -> Maybe (a, a)
support :: forall a. (Ord a, Num a) => Measure a -> Maybe (a, a)
support (Measure Piecewise (Poly a)
pieces) =
    case Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
pieces of
        [] -> Maybe (a, a)
forall a. Maybe a
Nothing
        [(Domain (Poly a), Poly a)]
ps -> (a, a) -> Maybe (a, a)
forall a. a -> Maybe a
Just ((a, Poly a) -> a
forall a b. (a, b) -> a
fst ((a, Poly a) -> a) -> (a, Poly a) -> a
forall a b. (a -> b) -> a -> b
$ [(a, Poly a)] -> (a, Poly a)
forall a. HasCallStack => [a] -> a
head [(a, Poly a)]
[(Domain (Poly a), Poly a)]
ps, (a, Poly a) -> a
forall a b. (a, b) -> a
fst ((a, Poly a) -> a) -> (a, Poly a) -> a
forall a b. (a -> b) -> a -> b
$ [(a, Poly a)] -> (a, Poly a)
forall a. HasCallStack => [a] -> a
last [(a, Poly a)]
[(Domain (Poly a), Poly a)]
ps)

-- | Check whether a signed measure is positive.
--
-- A signed measure is /positive/ if the measure of any set
-- is nonnegative. In other words a positive signed measure
-- is just a measure in the ordinary sense.
--
-- This test is nontrivial, as we have to check that the distribution
-- function is monotonically increasing.
isPositive :: (Ord a, Num a, Fractional a) => Measure a -> Bool
isPositive :: forall a. (Ord a, Num a, Fractional a) => Measure a -> Bool
isPositive (Measure Piecewise (Poly a)
m) = Poly a -> [(a, Poly a)] -> Bool
forall {a}.
(Ord a, Fractional a) =>
Poly a -> [(a, Poly a)] -> Bool
go Poly a
0 ([(a, Poly a)] -> Bool) -> [(a, Poly a)] -> Bool
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
m
  where
    go :: Poly a -> [(a, Poly a)] -> Bool
go Poly a
_ [] =
        Bool
True
    go Poly a
before ((a
x, Poly a
o) : []) =
        Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
before a
Domain (Poly a)
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
o a
Domain (Poly a)
x
    go Poly a
before ((a
x1, Poly a
o) : xos :: [(a, Poly a)]
xos@((a
x2, Poly a
_) : [(a, Poly a)]
_)) =
        (Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
before a
Domain (Poly a)
x1 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
o a
Domain (Poly a)
x1)
        Bool -> Bool -> Bool
&& Poly a -> (a, a) -> Bool
forall a. (Fractional a, Eq a, Ord a) => Poly a -> (a, a) -> Bool
Poly.isMonotonicallyIncreasingOn Poly a
o (a
x1,a
x2)
        Bool -> Bool -> Bool
&& Poly a -> [(a, Poly a)] -> Bool
go Poly a
o [(a, Poly a)]
xos

{-----------------------------------------------------------------------------
    Operations
    Numerical
------------------------------------------------------------------------------}
-- | Add two measures.
--
-- > total (add mx my) = total mx + total my
add :: (Ord a, Num a) => Measure a -> Measure a -> Measure a
add :: forall a. (Ord a, Num a) => Measure a -> Measure a -> Measure a
add (Measure Piecewise (Poly a)
mx) (Measure Piecewise (Poly a)
my) =
    Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a) -> Piecewise (Poly a)
forall a.
(Ord a, Num a) =>
Piecewise (Poly a) -> Piecewise (Poly a)
trim (Piecewise (Poly a) -> Piecewise (Poly a))
-> Piecewise (Poly a) -> Piecewise (Poly a)
forall a b. (a -> b) -> a -> b
$ (Poly a -> Poly a -> Poly a)
-> Piecewise (Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall o.
(Ord (Domain o), Num o) =>
(o -> o -> o) -> Piecewise o -> Piecewise o -> Piecewise o
Piecewise.zipPointwise Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
(+) Piecewise (Poly a)
mx Piecewise (Poly a)
my

-- | Scale a measure by a constant.
--
-- > total (scale a mx) = a * total mx
scale :: (Ord a, Num a) => a -> Measure a -> Measure a
scale :: forall a. (Ord a, Num a) => a -> Measure a -> Measure a
scale a
0 (Measure Piecewise (Poly a)
_) = Measure a
forall a. Num a => Measure a
zero
scale a
x (Measure Piecewise (Poly a)
m) = Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ (Poly a -> Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall o o'.
(Domain o ~ Domain o') =>
(o -> o') -> Piecewise o -> Piecewise o'
Piecewise.mapPieces (a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
Poly.scale a
x) Piecewise (Poly a)
m

-- | Translate a measure along the number line.
--
-- > eval (distribution (translate y m)) x
-- >    = eval (distribution m) (x - y)
translate :: (Ord a, Num a, Fractional a) => a -> Measure a -> Measure a
translate :: forall a.
(Ord a, Num a, Fractional a) =>
a -> Measure a -> Measure a
translate a
y (Measure Piecewise (Poly a)
m) =
    Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ (Domain (Poly a) -> Poly a -> Poly a)
-> Domain (Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall o.
(Ord (Domain o), Num (Domain o), Num o) =>
(Domain o -> o -> o) -> Domain o -> Piecewise o -> Piecewise o
Piecewise.translateWith a -> Poly a -> Poly a
Domain (Poly a) -> Poly a -> Poly a
forall a. (Fractional a, Eq a, Num a) => a -> Poly a -> Poly a
Poly.translate a
Domain (Poly a)
y Piecewise (Poly a)
m

{-----------------------------------------------------------------------------
    Operations
    Decomposition into continuous and discrete measures,
    needed for convolution.
------------------------------------------------------------------------------}
-- | Measure that is absolutely continuous
-- with respect to the Lebesgue measure,
-- Represented via its distribution function.
newtype Continuous a = Continuous { forall a. Continuous a -> Piecewise (Poly a)
unContinuous :: Piecewise (Poly a) }
    -- INVARIANT: The last piece is @Poly.constant p@ for some @p :: a@.

-- | Density function (Radon–Nikodym derivative) of an absolutely
-- continuous measure.
newtype Density a = Density (Piecewise (Poly a))
    -- INVARIANT: The last piece is @Poly.constant 0@.

-- | Density function of an absolutely continuous measure.
toDensity
    :: (Ord a, Num a, Fractional a)
    => Continuous a -> Density a
toDensity :: forall a. (Ord a, Num a, Fractional a) => Continuous a -> Density a
toDensity = Piecewise (Poly a) -> Density a
forall a. Piecewise (Poly a) -> Density a
Density (Piecewise (Poly a) -> Density a)
-> (Continuous a -> Piecewise (Poly a))
-> Continuous a
-> Density a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Poly a -> Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall o o'.
(Domain o ~ Domain o') =>
(o -> o') -> Piecewise o -> Piecewise o'
Piecewise.mapPieces Poly a -> Poly a
forall a. Num a => Poly a -> Poly a
Poly.differentiate (Piecewise (Poly a) -> Piecewise (Poly a))
-> (Continuous a -> Piecewise (Poly a))
-> Continuous a
-> Piecewise (Poly a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Continuous a -> Piecewise (Poly a)
forall a. Continuous a -> Piecewise (Poly a)
unContinuous

-- | Decompose a mixed measure into
-- a continuous measure and a discrete measure.
-- See also [Lebesgue's decomposition theorem
-- ](https://en.wikipedia.org/wiki/Lebesgue%27s_decomposition_theorem)
decompose
    :: (Ord a, Num a, Fractional a)
    => Measure a -> (Continuous a, D.Discrete a)
decompose :: forall a.
(Ord a, Num a, Fractional a) =>
Measure a -> (Continuous a, Discrete a)
decompose (Measure Piecewise (Poly a)
m) =
    ( Piecewise (Poly a) -> Continuous a
forall a. Piecewise (Poly a) -> Continuous a
Continuous (Piecewise (Poly a) -> Continuous a)
-> Piecewise (Poly a) -> Continuous a
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a) -> Piecewise (Poly a)
forall a.
(Ord a, Num a) =>
Piecewise (Poly a) -> Piecewise (Poly a)
trim (Piecewise (Poly a) -> Piecewise (Poly a))
-> Piecewise (Poly a) -> Piecewise (Poly a)
forall a b. (a -> b) -> a -> b
$ [(Domain (Poly a), Poly a)] -> Piecewise (Poly a)
forall o. Ord (Domain o) => [(Domain o, o)] -> Piecewise o
Piecewise.fromAscPieces [(a, Poly a)]
[(Domain (Poly a), Poly a)]
withoutJumps
    , Map a a -> Discrete a
forall a. (Ord a, Num a) => Map a a -> Discrete a
D.fromMap (Map a a -> Discrete a) -> Map a a -> Discrete a
forall a b. (a -> b) -> a -> b
$ [(a, a)] -> Map a a
forall k a. Ord k => [(k, a)] -> Map k a
Map.fromList [(a, a)]
jumps
    )
  where
    pieces :: [(Domain (Poly a), Poly a)]
pieces = Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
m

    withoutJumps :: [(a, Poly a)]
withoutJumps =
        ((a, Poly a) -> a -> (a, Poly a))
-> [(a, Poly a)] -> [a] -> [(a, Poly a)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\(a
x,Poly a
o) a
j -> (a
x, Poly a
o Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
- a -> Poly a
forall a. a -> Poly a
Poly.constant a
j)) [(a, Poly a)]
pieces [a]
totalJumps
    totalJumps :: [a]
totalJumps = [a] -> [a]
forall a. HasCallStack => [a] -> [a]
tail ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ (a -> a -> a) -> a -> [a] -> [a]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl' a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0 ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ ((a, a) -> a) -> [(a, a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a, a) -> a
forall a b. (a, b) -> b
snd [(a, a)]
jumps

    jumps :: [(a, a)]
jumps = Poly a -> [(a, Poly a)] -> [(a, a)]
forall {b}. Num b => Poly b -> [(b, Poly b)] -> [(b, b)]
go Poly a
0 [(a, Poly a)]
pieces
      where
        go :: Poly b -> [(b, Poly b)] -> [(b, b)]
go Poly b
_ [] = []
        go Poly b
prev ((b
x,Poly b
o) : [(b, Poly b)]
xos) =
            (b
x, Poly b -> b -> b
forall a. Num a => Poly a -> a -> a
Poly.eval Poly b
o b
x b -> b -> b
forall a. Num a => a -> a -> a
- Poly b -> b -> b
forall a. Num a => Poly a -> a -> a
Poly.eval Poly b
prev b
x) (b, b) -> [(b, b)] -> [(b, b)]
forall a. a -> [a] -> [a]
: Poly b -> [(b, Poly b)] -> [(b, b)]
go Poly b
o [(b, Poly b)]
xos

{-----------------------------------------------------------------------------
    Observations
    Integration
------------------------------------------------------------------------------}
-- | Integrate a polynomial @f@ with respect to the given measure @m@,
-- \( \int f(x) dm(x) \).
integrate :: (Ord a, Num a, Fractional a) => Poly a -> Measure a -> a
integrate :: forall a. (Ord a, Num a, Fractional a) => Poly a -> Measure a -> a
integrate Poly a
f Measure a
m =
    Poly a -> Continuous a -> a
forall a.
(Ord a, Num a, Fractional a) =>
Poly a -> Continuous a -> a
integrateContinuous Poly a
f Continuous a
continuous
    a -> a -> a
forall a. Num a => a -> a -> a
+ (a -> a) -> Discrete a -> a
forall a. (Ord a, Num a) => (a -> a) -> Discrete a -> a
D.integrate (Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
f) Discrete a
discrete
  where
    (Continuous a
continuous, Discrete a
discrete) = Measure a -> (Continuous a, Discrete a)
forall a.
(Ord a, Num a, Fractional a) =>
Measure a -> (Continuous a, Discrete a)
decompose Measure a
m

-- | Integrate a polynomial over an absolutely continuous measure.
integrateContinuous
    :: (Ord a, Num a, Fractional a)
    => Poly a -> Continuous a -> a
integrateContinuous :: forall a.
(Ord a, Num a, Fractional a) =>
Poly a -> Continuous a -> a
integrateContinuous Poly a
f Continuous a
gg
    | [(a, Poly a)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, Poly a)]
gpieces = a
0
    | Bool
otherwise = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ (((a, a), Poly a) -> a) -> [((a, a), Poly a)] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ((a, a), Poly a) -> a
forall {a}. (Eq a, Fractional a) => ((a, a), Poly a) -> a
integrateOverInterval ([((a, a), Poly a)] -> [a]) -> [((a, a), Poly a)] -> [a]
forall a b. (a -> b) -> a -> b
$ [((a, a), Poly a)]
integrands
  where
    Density Piecewise (Poly a)
g = Continuous a -> Density a
forall a. (Ord a, Num a, Fractional a) => Continuous a -> Density a
toDensity Continuous a
gg
    gpieces :: [(Domain (Poly a), Poly a)]
gpieces = Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
g

    -- Pieces on the bounded intervals
    boundedPieces :: [(b, b)] -> [((b, b), b)]
boundedPieces [(b, b)]
xos =
        ((b, b) -> (b, b) -> ((b, b), b))
-> [(b, b)] -> [(b, b)] -> [((b, b), b)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\(b
x1,b
o) (b
x2,b
_) -> ((b
x1, b
x2), b
o)) [(b, b)]
xos (Int -> [(b, b)] -> [(b, b)]
forall a. Int -> [a] -> [a]
drop Int
1 [(b, b)]
xos)

    integrands :: [((a, a), Poly a)]
integrands = [ ((a, a)
x12, Poly a
f Poly a -> Poly a -> Poly a
forall a. Num a => a -> a -> a
* Poly a
o) | ((a, a)
x12, Poly a
o) <- [(a, Poly a)] -> [((a, a), Poly a)]
forall {b} {b}. [(b, b)] -> [((b, b), b)]
boundedPieces [(a, Poly a)]
gpieces ]

    integrateOverInterval :: ((a, a), Poly a) -> a
integrateOverInterval ((a
x1, a
x2), Poly a
p) =
        Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
pp a
Domain (Poly a)
x2 a -> a -> a
forall a. Num a => a -> a -> a
- Poly a -> Domain (Poly a) -> Codomain (Poly a)
forall f. Function f => f -> Domain f -> Codomain f
eval Poly a
pp a
Domain (Poly a)
x1
      where
        pp :: Poly a
pp = Poly a -> Poly a
forall a. (Eq a, Fractional a) => Poly a -> Poly a
Poly.integrate Poly a
p

{-----------------------------------------------------------------------------
    Operations
    Convolution
------------------------------------------------------------------------------}
{-$ NOTE [Convolution]

In order to compute a convolution,
we convolve a density with the distribution function.

Let $f$ denote a density, which can be continuous or a Dirac delta.
Let $G$ denote a distribution function.
Let $H = f * G$ be the result of the convolution.
It can be shown that this is the distribution function of the
convolution of the densities, $h = f * g$.

The formula for convolution is

$ H(y) = ∫ f(y - x) G(x) dx = ∫ f (x) G(y - x) dx$.

When $f$ is a sum of delta functions, $f = Σ w_j delta_{x_j}(x)$,
this integral becomes ($y - x = x_j$ => $x = y - x_j$)

$ H(y) = Σ w_j G(y - x_j) $.

When $f$ is a piecewise polynomial, we can convolve the pieces.

When convolving with a distribution function, the final piece
will be a constant $g_n$ on the interval $[x_n,∞)$.
In this case, the convolution is given by

\[
H(y)
    = ∫ f (x) G(y - x) dx
    = ∫_{ -∞}^{y-x_n} f(x) g_n dx
    = g_n F(y-x_n)
\]

where $F$ is the distribution function of the density $f$.
-}

-- | Convolve a discrete measure with a mixed measure.
--
-- See NOTE [Convolution].
convolveDiscrete
    :: (Ord a, Num a, Fractional a)
    => D.Discrete a -> Measure a -> Measure a
convolveDiscrete :: forall a.
(Ord a, Num a, Fractional a) =>
Discrete a -> Measure a -> Measure a
convolveDiscrete Discrete a
f Measure a
gg =
    (Measure a -> Measure a -> Measure a)
-> Measure a -> [Measure a] -> Measure a
forall a b. (a -> b -> b) -> b -> [a] -> b
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr Measure a -> Measure a -> Measure a
forall a. (Ord a, Num a) => Measure a -> Measure a -> Measure a
add Measure a
forall a. Num a => Measure a
zero
        [ a -> Measure a -> Measure a
forall a. (Ord a, Num a) => a -> Measure a -> Measure a
scale a
w (a -> Measure a -> Measure a
forall a.
(Ord a, Num a, Fractional a) =>
a -> Measure a -> Measure a
translate a
x Measure a
gg)
        | (a
x, a
w) <- Map a a -> [(a, a)]
forall k a. Map k a -> [(k, a)]
Map.toAscList (Map a a -> [(a, a)]) -> Map a a -> [(a, a)]
forall a b. (a -> b) -> a -> b
$ Discrete a -> Map a a
forall a. Num a => Discrete a -> Map a a
D.toMap Discrete a
f
        ]

-- | Convolve an absolutely continuous measure with a mixed measure.
--
-- See NOTE [Convolution].
convolveContinuous
    :: (Ord a, Num a, Fractional a)
    => Continuous a -> Measure a -> Measure a
convolveContinuous :: forall a.
(Ord a, Num a, Fractional a) =>
Continuous a -> Measure a -> Measure a
convolveContinuous (Continuous Piecewise (Poly a)
ff) (Measure Piecewise (Poly a)
gg)
    | [(a, Poly a)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, Poly a)]
ffpieces = Measure a
forall a. Num a => Measure a
zero
    | [(a, Poly a)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(a, Poly a)]
ggpieces = Measure a
forall a. Num a => Measure a
zero
    | Bool
otherwise = Piecewise (Poly a) -> Measure a
forall a. Piecewise (Poly a) -> Measure a
Measure (Piecewise (Poly a) -> Measure a)
-> Piecewise (Poly a) -> Measure a
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a) -> Piecewise (Poly a)
forall a.
(Ord a, Num a) =>
Piecewise (Poly a) -> Piecewise (Poly a)
trim (Piecewise (Poly a) -> Piecewise (Poly a))
-> Piecewise (Poly a) -> Piecewise (Poly a)
forall a b. (a -> b) -> a -> b
$ Piecewise (Poly a)
boundedConvolutions Piecewise (Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall a. Num a => a -> a -> a
+ Piecewise (Poly a)
lastConvolution
  where
    ffpieces :: [(Domain (Poly a), Poly a)]
ffpieces = Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
ff
    ggpieces :: [(Domain (Poly a), Poly a)]
ggpieces = Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
gg

    Density Piecewise (Poly a)
f = Continuous a -> Density a
forall a. (Ord a, Num a, Fractional a) => Continuous a -> Density a
toDensity (Piecewise (Poly a) -> Continuous a
forall a. Piecewise (Poly a) -> Continuous a
Continuous Piecewise (Poly a)
ff)
    fpieces :: [(Domain (Poly a), Poly a)]
fpieces = Piecewise (Poly a) -> [(Domain (Poly a), Poly a)]
forall o. Ord (Domain o) => Piecewise o -> [(Domain o, o)]
Piecewise.toAscPieces Piecewise (Poly a)
f

    -- Pieces on the bounded intervals
    boundedPieces :: [(b, b)] -> [(b, b, b)]
boundedPieces [(b, b)]
xos =
        ((b, b) -> (b, b) -> (b, b, b))
-> [(b, b)] -> [(b, b)] -> [(b, b, b)]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\(b
x,b
o) (b
y,b
_) -> (b
x, b
y, b
o)) [(b, b)]
xos (Int -> [(b, b)] -> [(b, b)]
forall a. Int -> [a] -> [a]
drop Int
1 [(b, b)]
xos)

    boundedConvolutions :: Piecewise (Poly a)
boundedConvolutions =
        [Piecewise (Poly a)] -> Piecewise (Poly a)
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Piecewise (Poly a)] -> Piecewise (Poly a))
-> [Piecewise (Poly a)] -> Piecewise (Poly a)
forall a b. (a -> b) -> a -> b
$
            [ [(Domain (Poly a), Poly a)] -> Piecewise (Poly a)
forall o. Ord (Domain o) => [(Domain o, o)] -> Piecewise o
Piecewise.fromAscPieces ((a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]
forall a.
(Fractional a, Eq a, Ord a) =>
(a, a, Poly a) -> (a, a, Poly a) -> [(a, Poly a)]
Poly.convolve (a, a, Poly a)
fo (a, a, Poly a)
ggo)
            | (a, a, Poly a)
fo <- [(a, Poly a)] -> [(a, a, Poly a)]
forall {b} {b}. [(b, b)] -> [(b, b, b)]
boundedPieces [(a, Poly a)]
fpieces
            , (a, a, Poly a)
ggo <- [(a, Poly a)] -> [(a, a, Poly a)]
forall {b} {b}. [(b, b)] -> [(b, b, b)]
boundedPieces [(a, Poly a)]
ggpieces
            ]

    (a
xlast, Poly a
plast) = [(a, Poly a)] -> (a, Poly a)
forall a. HasCallStack => [a] -> a
last [(a, Poly a)]
ggpieces
    glast :: a
glast = case Poly a -> [a]
forall a. Poly a -> [a]
Poly.toCoefficients Poly a
plast of
        [] -> a
0
        (a
a0:[a]
_) -> a
a0
    lastConvolution :: Piecewise (Poly a)
lastConvolution =
        (Poly a -> Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall o o'.
(Domain o ~ Domain o') =>
(o -> o') -> Piecewise o -> Piecewise o'
Piecewise.mapPieces (a -> Poly a -> Poly a
forall a. Num a => a -> Poly a -> Poly a
Poly.scale a
glast)
        (Piecewise (Poly a) -> Piecewise (Poly a))
-> Piecewise (Poly a) -> Piecewise (Poly a)
forall a b. (a -> b) -> a -> b
$ (Domain (Poly a) -> Poly a -> Poly a)
-> Domain (Poly a) -> Piecewise (Poly a) -> Piecewise (Poly a)
forall o.
(Ord (Domain o), Num (Domain o), Num o) =>
(Domain o -> o -> o) -> Domain o -> Piecewise o -> Piecewise o
Piecewise.translateWith a -> Poly a -> Poly a
Domain (Poly a) -> Poly a -> Poly a
forall a. (Fractional a, Eq a, Num a) => a -> Poly a -> Poly a
Poly.translate a
Domain (Poly a)
xlast Piecewise (Poly a)
ff

-- | Additive convolution of two measures.
--
-- Properties:
--
-- > convolve (dirac x) (dirac y) = dirac (x + y)
-- >
-- > convolve mx my               =  convolve my mx
-- > convolve (add mx my) mz      =  add (convolve mx mz) (convolve my mz)
-- > translate z (convolve mx my) =  convolve (translate z mx) my
-- > total (convolve mx my)       =  total mx * total my
convolve
    :: (Ord a, Num a, Fractional a)
    => Measure a -> Measure a -> Measure a
convolve :: forall a.
(Ord a, Num a, Fractional a) =>
Measure a -> Measure a -> Measure a
convolve Measure a
mx Measure a
my =
    Measure a -> Measure a -> Measure a
forall a. (Ord a, Num a) => Measure a -> Measure a -> Measure a
add (Continuous a -> Measure a -> Measure a
forall a.
(Ord a, Num a, Fractional a) =>
Continuous a -> Measure a -> Measure a
convolveContinuous Continuous a
contx Measure a
my) (Discrete a -> Measure a -> Measure a
forall a.
(Ord a, Num a, Fractional a) =>
Discrete a -> Measure a -> Measure a
convolveDiscrete Discrete a
deltax Measure a
my)
  where
    (Continuous a
contx, Discrete a
deltax) = Measure a -> (Continuous a, Discrete a)
forall a.
(Ord a, Num a, Fractional a) =>
Measure a -> (Continuous a, Discrete a)
decompose Measure a
mx