{-# LANGUAGE CPP #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE NoImplicitPrelude #-}

-- |
-- Module      : Numeric.Uncertain.MonteCarlo
-- Copyright   : (c) Justin Le 2016
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- Provides an interface for computing and propagating uncertainty by using
-- <https://en.wikipedia.org/wiki/Monte_Carlo_method Monte Carlo simulations>.
-- These may be more accurate than the normal propagators if the third-order
-- taylor series expansion terms are non-negligible.
--
-- Basically simulates sampling from the distribution represented by the given
-- 'Uncert's, applying the function of interest, and aggregating the mean
-- and standard deviation of the results.  @x '+/-' dx@ is treated as
-- a random variable whose probability density is the normal distribution
-- with mean @x@ and standard deviation @dx@.
--
-- This module attempts to duplicate the API offered by "Numeric.Uncertain"
-- and is meant to be imported qualified alongside "Numeric.Uncertain"
--
-- @
-- import           Numeric.Uncertain
-- import qualified Numeric.Uncertain.MonteCarlo as MC
-- @
--
-- Actions are parameterized over all 'PrimMonad' instances, so can be run
-- under both 'ST' and 'IO', making it suitable for exploratory purposes.
-- All functions require a 'Gen' from "System.Random.MWC" for random value
-- generation purposes.
--
-- @
-- ghci> import qualified Numeric.Uncertain.MonteCarlo as MC
-- ghci> import System.Random.MWC
-- ghci> let x = 1.52 '+/-' 0.07
-- ghci> let y = 781.4 +/- 0.3
-- ghci> let z = 1.53e-1 `'withPrecision'` 3
-- ghci> g <- 'create'
-- ghci> cosh x
-- 2.4 +/- 0.2
-- ghci> MC.liftU cosh x g
-- 2.4 +/- 0.2
-- ghci> exp x / z * sin (y ** z)
-- 10.9 +/- 0.9
-- ghci> MC.liftU3 (\\a b c -> exp a / c * sin (b**c)) x y z g
-- 10.8 +/- 1.0
-- ghci> pi + 3 * logBase x y
-- 52 +/- 5
-- ghci> MC.liftU2 (\\a b -> pi + 3 * logBase a b) x y g
-- 51 +/- 5
-- @
module Numeric.Uncertain.MonteCarlo (
  -- * Sampling from an 'Uncert'
  sampleUncert,

  -- * Lifting functions via Monte Carlo simulation

  -- ** Fixed iterations
  liftU,
  liftU2,
  liftU3,
  liftU4,
  liftU5,
  liftUF,

  -- ** Variable iterations
  liftU',
  liftU2',
  liftU3',
  liftU4',
  liftU5',
  liftUF',
)
where

import Control.Monad
import Control.Monad.Primitive
import Data.Hople
import Numeric.Uncertain (Uncert, fromSamples, uMeanStd)
import Prelude.Compat
import System.Random.MWC
import System.Random.MWC.Distributions

-- | Sample a random 'Double' from the distribution specified by an
-- @'Uncert' 'Double'@.  @x '+/-' dx@ is treated as a random variable whose
-- probability density is the normal distribution with mean @x@ and
-- standard deviation @dx@.
sampleUncert ::
  PrimMonad m =>
  Uncert Double ->
  Gen (PrimState m) ->
  m Double
sampleUncert :: forall (m :: * -> *).
PrimMonad m =>
Uncert Double -> Gen (PrimState m) -> m Double
sampleUncert (Uncert Double -> (Double, Double)
forall a. Floating a => Uncert a -> (a, a)
uMeanStd -> (Double
x, Double
dx)) = Double -> Double -> Gen (PrimState m) -> m Double
forall g (m :: * -> *).
StatefulGen g m =>
Double -> Double -> g -> m Double
normal Double
x Double
dx
{-# INLINEABLE sampleUncert #-}

-- | Lifts a numeric function over an 'Uncert' using a Monte Carlo
-- simulation with 1000 samples.
--
-- @
-- ghci> g <- 'create'
-- ghci> MC.liftU (\\x -> log x ^ 2) (12.2 +/- 0.5) g
-- 6.3 +/- 0.2
-- @
liftU ::
  PrimMonad m =>
  (Double -> Double) ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU :: forall (m :: * -> *).
PrimMonad m =>
(Double -> Double)
-> Uncert Double -> Gen (PrimState m) -> m (Uncert Double)
liftU = Int
-> (Double -> Double)
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double)
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU' Int
1000
{-# INLINE liftU #-}

-- | Lifts a multivariate numeric function on a container (given as an @f
-- a -> a@) to work on a container of 'Uncert's using a Monte Carlo
-- simulation with 1000 samples.
--
-- @
-- ghci> g <- 'create'
-- ghci> M.liftUF (\\[x,y,z] -> x*y+z) [12.2 +/- 0.5, 56 +/- 2, 0.12 +/- 0.08] g
-- 680 +/- 40
-- @
liftUF ::
  (Traversable f, PrimMonad m) =>
  (f Double -> Double) ->
  f (Uncert Double) ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftUF :: forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
(f Double -> Double)
-> f (Uncert Double) -> Gen (PrimState m) -> m (Uncert Double)
liftUF = Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
liftUF' Int
1000
{-# INLINE liftUF #-}

-- | Lifts a two-argument (curried) function over two 'Uncert's using
-- a Monte Carlo simulation with 1000 samples.
--
-- @
-- ghci> g <- 'create'
-- ghci> MC.liftU2 (\\x y -> x**y) (13.5 +/- 0.1) (1.64 +/- 0.08)
-- 70 +/- 20
-- @
liftU2 ::
  PrimMonad m =>
  (Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU2 :: forall (m :: * -> *).
PrimMonad m =>
(Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU2 = Int
-> (Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU2' Int
1000
{-# INLINE liftU2 #-}

-- | Lifts a three-argument (curried) function over three 'Uncert's.  See
-- 'liftU2' and 'liftUF' for more details.
liftU3 ::
  PrimMonad m =>
  (Double -> Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU3 :: forall (m :: * -> *).
PrimMonad m =>
(Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU3 = Int
-> (Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU3' Int
1000
{-# INLINE liftU3 #-}

-- | Lifts a four-argument (curried) function over four 'Uncert's.  See
-- 'liftU2' and 'liftUF' for more details.
liftU4 ::
  (PrimMonad m, Applicative m) =>
  (Double -> Double -> Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU4 :: forall (m :: * -> *).
(PrimMonad m, Applicative m) =>
(Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU4 = Int
-> (Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU4' Int
1000
{-# INLINE liftU4 #-}

-- | Lifts a five-argument (curried) function over five 'Uncert's.  See
-- 'liftU2' and 'liftUF' for more details.
liftU5 ::
  PrimMonad m =>
  (Double -> Double -> Double -> Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU5 :: forall (m :: * -> *).
PrimMonad m =>
(Double -> Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU5 = Int
-> (Double -> Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU5' Int
1000
{-# INLINE liftU5 #-}

-- | Like 'liftU', but allows you to specify the number of samples to run
-- the Monte Carlo simulation with.
liftU' ::
  PrimMonad m =>
  Int ->
  (Double -> Double) ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU' :: forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double)
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU' Int
n Double -> Double
f Uncert Double
u Gen (PrimState m)
g = [Double] -> Uncert Double
forall a. Fractional a => [a] -> Uncert a
fromSamples ([Double] -> Uncert Double) -> m [Double] -> m (Uncert Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> m Double -> m [Double]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM Int
n m Double
samp
  where
    samp :: m Double
samp = Double -> Double
f (Double -> Double) -> m Double -> m Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Uncert Double -> Gen (PrimState m) -> m Double
forall (m :: * -> *).
PrimMonad m =>
Uncert Double -> Gen (PrimState m) -> m Double
sampleUncert Uncert Double
u Gen (PrimState m)
g
{-# INLINEABLE liftU' #-}

-- | Like 'liftUF', but allows you to specify the number of samples to run
-- the Monte Carlo simulation with.
liftUF' ::
  (Traversable f, PrimMonad m) =>
  Int ->
  (f Double -> Double) ->
  f (Uncert Double) ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftUF' :: forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
liftUF' Int
n f Double -> Double
f f (Uncert Double)
us Gen (PrimState m)
g = [Double] -> Uncert Double
forall a. Fractional a => [a] -> Uncert a
fromSamples ([Double] -> Uncert Double) -> m [Double] -> m (Uncert Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> m Double -> m [Double]
forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM Int
n m Double
samp
  where
    samp :: m Double
samp = f Double -> Double
f (f Double -> Double) -> m (f Double) -> m Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Uncert Double -> m Double) -> f (Uncert Double) -> m (f Double)
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> f a -> f (f b)
traverse (Uncert Double -> Gen (PrimState m) -> m Double
forall (m :: * -> *).
PrimMonad m =>
Uncert Double -> Gen (PrimState m) -> m Double
`sampleUncert` Gen (PrimState m)
g) f (Uncert Double)
us
{-# INLINEABLE liftUF' #-}

-- | Like 'liftU2', but allows you to specify the number of samples to run
-- the Monte Carlo simulation with.
liftU2' ::
  PrimMonad m =>
  Int ->
  (Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU2' :: forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU2' Int
n Double -> Double -> Double
f Uncert Double
x Uncert Double
y = Int
-> (H2 Double -> Double)
-> H2 (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
liftUF' Int
n ((Double -> Double -> Double) -> H2 Double -> Double
forall a. (a -> a -> a) -> H2 a -> a
uncurryH2 Double -> Double -> Double
f) (Uncert Double -> Uncert Double -> H2 (Uncert Double)
forall a. a -> a -> H2 a
H2 Uncert Double
x Uncert Double
y)
{-# INLINEABLE liftU2' #-}

-- | Like 'liftU3', but allows you to specify the number of samples to run
-- the Monte Carlo simulation with.
liftU3' ::
  PrimMonad m =>
  Int ->
  (Double -> Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU3' :: forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU3' Int
n Double -> Double -> Double -> Double
f Uncert Double
x Uncert Double
y Uncert Double
z = Int
-> (H3 Double -> Double)
-> H3 (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
liftUF' Int
n ((Double -> Double -> Double -> Double) -> H3 Double -> Double
forall a. (a -> a -> a -> a) -> H3 a -> a
uncurryH3 Double -> Double -> Double -> Double
f) (Uncert Double
-> Uncert Double -> Uncert Double -> H3 (Uncert Double)
forall a. a -> a -> a -> H3 a
H3 Uncert Double
x Uncert Double
y Uncert Double
z)
{-# INLINEABLE liftU3' #-}

-- | Like 'liftU4', but allows you to specify the number of samples to run
-- the Monte Carlo simulation with.
liftU4' ::
  PrimMonad m =>
  Int ->
  (Double -> Double -> Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU4' :: forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU4' Int
n Double -> Double -> Double -> Double -> Double
f Uncert Double
x Uncert Double
y Uncert Double
z Uncert Double
a = Int
-> (H4 Double -> Double)
-> H4 (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
liftUF' Int
n ((Double -> Double -> Double -> Double -> Double)
-> H4 Double -> Double
forall a. (a -> a -> a -> a -> a) -> H4 a -> a
uncurryH4 Double -> Double -> Double -> Double -> Double
f) (Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> H4 (Uncert Double)
forall a. a -> a -> a -> a -> H4 a
H4 Uncert Double
x Uncert Double
y Uncert Double
z Uncert Double
a)
{-# INLINEABLE liftU4' #-}

-- | Like 'liftU5', but allows you to specify the number of samples to run
-- the Monte Carlo simulation with.
liftU5' ::
  PrimMonad m =>
  Int ->
  (Double -> Double -> Double -> Double -> Double -> Double) ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Uncert Double ->
  Gen (PrimState m) ->
  m (Uncert Double)
liftU5' :: forall (m :: * -> *).
PrimMonad m =>
Int
-> (Double -> Double -> Double -> Double -> Double -> Double)
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Gen (PrimState m)
-> m (Uncert Double)
liftU5' Int
n Double -> Double -> Double -> Double -> Double -> Double
f Uncert Double
x Uncert Double
y Uncert Double
z Uncert Double
a Uncert Double
b = Int
-> (H5 Double -> Double)
-> H5 (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
forall (f :: * -> *) (m :: * -> *).
(Traversable f, PrimMonad m) =>
Int
-> (f Double -> Double)
-> f (Uncert Double)
-> Gen (PrimState m)
-> m (Uncert Double)
liftUF' Int
n ((Double -> Double -> Double -> Double -> Double -> Double)
-> H5 Double -> Double
forall a. (a -> a -> a -> a -> a -> a) -> H5 a -> a
uncurryH5 Double -> Double -> Double -> Double -> Double -> Double
f) (Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> Uncert Double
-> H5 (Uncert Double)
forall a. a -> a -> a -> a -> a -> H5 a
H5 Uncert Double
x Uncert Double
y Uncert Double
z Uncert Double
a Uncert Double
b)
{-# INLINEABLE liftU5' #-}