-- | Random useful stuff I didn't know where to put.
module Bio.Util.Numeric (
    wilson, invnormcdf, choose,
    estimateComplexity, showNum, showOOM, readOOM,
    log1p, expm1, (<#>),
    log1mexp, log1pexp,
    lsum, llerp
                ) where

import Prelude
import Data.Char ( intToDigit )
import Data.List ( foldl1' )
import Options.Applicative.Builder ( eitherReader, ReadM )

-- | Calculates the Wilson Score interval.
-- If @(l,m,h) = wilson c x n@, then @m@ is the binary proportion and
-- @(l,h)@ it's @c@-confidence interval for @x@ positive examples out of
-- @n@ observations.  @c@ is typically something like 0.05.

wilson :: Double -> Int -> Int -> (Double, Double, Double)
wilson :: Double -> Int -> Int -> (Double, Double, Double)
wilson c :: Double
c x :: Int
x n :: Int
n = ( (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
h) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d, Double
p, (Double
m Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
d )
  where
    nn :: Double
nn = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
    p :: Double
p  = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nn

    z :: Double
z = Double -> Double
forall a. (Ord a, Floating a) => a -> a
invnormcdf (1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*0.5)
    h :: Double
h = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (( Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* (1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 0.25Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nn ) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nn)
    m :: Double
m = Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nn
    d :: Double
d = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
nn

showNum :: Show a => a -> String
showNum :: a -> String
showNum = String -> String -> String
triplets [] (String -> String) -> (a -> String) -> a -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> String
forall a. [a] -> [a]
reverse (String -> String) -> (a -> String) -> a -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> String
forall a. Show a => a -> String
show
  where
    triplets :: String -> String -> String
triplets acc :: String
acc [] = String
acc
    triplets acc :: String
acc [a :: Char
a] = Char
aChar -> String -> String
forall a. a -> [a] -> [a]
:String
acc
    triplets acc :: String
acc [a :: Char
a,b :: Char
b] = Char
bChar -> String -> String
forall a. a -> [a] -> [a]
:Char
aChar -> String -> String
forall a. a -> [a] -> [a]
:String
acc
    triplets acc :: String
acc [a :: Char
a,b :: Char
b,c :: Char
c] = Char
cChar -> String -> String
forall a. a -> [a] -> [a]
:Char
bChar -> String -> String
forall a. a -> [a] -> [a]
:Char
aChar -> String -> String
forall a. a -> [a] -> [a]
:String
acc
    triplets acc :: String
acc (a :: Char
a:b :: Char
b:c :: Char
c:s :: String
s) = String -> String -> String
triplets (','Char -> String -> String
forall a. a -> [a] -> [a]
:Char
cChar -> String -> String
forall a. a -> [a] -> [a]
:Char
bChar -> String -> String
forall a. a -> [a] -> [a]
:Char
aChar -> String -> String
forall a. a -> [a] -> [a]
:String
acc) String
s


showOOM :: (Enum a, Num a, Ord a) => a -> String
showOOM :: a -> String
showOOM x :: a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0     = '-' Char -> String -> String
forall a. a -> [a] -> [a]
: a -> String
forall a. (Enum a, Num a, Ord a) => a -> String
showOOM (a -> a
forall a. Num a => a -> a
negate a
x)
          | Bool
otherwise = Int -> String -> String
findSuffix (a -> Int
forall a. Enum a => a -> Int
fromEnum (a
x a -> a -> a
forall a. Num a => a -> a -> a
* 100 a -> a -> a
forall a. Num a => a -> a -> a
+ 5) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 10) ".kMGTPEZY"
  where
    findSuffix :: Int -> String -> String
findSuffix _ [    ] = "many"
    findSuffix y :: Int
y (s :: Char
s:ss :: String
ss) | Int
y Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 100  = Int -> Char
intToDigit (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
y 10) Char -> String -> String
forall a. a -> [a] -> [a]
: case (Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
y 10, Char
s) of
                                            (0,'.') -> [] ; (0,_) -> [Char
s] ; (d :: Int
d,_) -> [Char
s, Int -> Char
intToDigit Int
d]
                        | Int
y Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 1000 = Int -> Char
intToDigit (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
y 100) Char -> String -> String
forall a. a -> [a] -> [a]
: Int -> Char
intToDigit (Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
y 100 Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 10) Char -> String -> String
forall a. a -> [a] -> [a]
:
                                            if Char
s Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '.' then [] else [Char
s]
                        | Int
y Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 10000 = Int -> Char
intToDigit (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
y 1000) Char -> String -> String
forall a. a -> [a] -> [a]
: Int -> Char
intToDigit (Int -> Int -> Int
forall a. Integral a => a -> a -> a
mod Int
y 1000 Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 100) Char -> String -> String
forall a. a -> [a] -> [a]
:
                                            '0' Char -> String -> String
forall a. a -> [a] -> [a]
: if Char
s Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '.' then [] else [Char
s]
                        | Bool
otherwise = Int -> String -> String
findSuffix (Int -> Int -> Int
forall a. Integral a => a -> a -> a
div Int
y 1000) String
ss

readOOM :: (Read a, Num a) => ReadM a
readOOM :: ReadM a
readOOM = (String -> Either String a) -> ReadM a
forall a. (String -> Either String a) -> ReadM a
eitherReader ((String -> Either String a) -> ReadM a)
-> (String -> Either String a) -> ReadM a
forall a b. (a -> b) -> a -> b
$ \s :: String
s -> case ReadS a
forall a. Read a => ReadS a
reads String
s of
    [(n :: a
n,[ ])] -> a -> Either String a
forall a b. b -> Either a b
Right a
n
    [(n :: a
n,"k")] -> a -> Either String a
forall a b. b -> Either a b
Right (a -> Either String a) -> a -> Either String a
forall a b. (a -> b) -> a -> b
$ a
n a -> a -> a
forall a. Num a => a -> a -> a
* 1000
    [(n :: a
n,"M")] -> a -> Either String a
forall a b. b -> Either a b
Right (a -> Either String a) -> a -> Either String a
forall a b. (a -> b) -> a -> b
$ a
n a -> a -> a
forall a. Num a => a -> a -> a
* 1000000
    [(n :: a
n,"G")] -> a -> Either String a
forall a b. b -> Either a b
Right (a -> Either String a) -> a -> Either String a
forall a b. (a -> b) -> a -> b
$ a
n a -> a -> a
forall a. Num a => a -> a -> a
* 1000000000
    [(n :: a
n,"T")] -> a -> Either String a
forall a b. b -> Either a b
Right (a -> Either String a) -> a -> Either String a
forall a b. (a -> b) -> a -> b
$ a
n a -> a -> a
forall a. Num a => a -> a -> a
* 1000000000000
    _         -> String -> Either String a
forall a b. a -> Either a b
Left (String -> Either String a) -> String -> Either String a
forall a b. (a -> b) -> a -> b
$ "unable to parse: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ String -> String
forall a. Show a => a -> String
show String
s


-- Stolen from Lennart Augustsson's erf package, who in turn took it from
-- <http://home.online.no/~pjacklam/notes/invnorm/> Accurate to about 1e-9.
invnormcdf :: (Ord a, Floating a) => a -> a
invnormcdf :: a -> a
invnormcdf p :: a
p =
    let a1 :: a
a1 = -3.969683028665376e+01
        a2 :: a
a2 =  2.209460984245205e+02
        a3 :: a
a3 = -2.759285104469687e+02
        a4 :: a
a4 =  1.383577518672690e+02
        a5 :: a
a5 = -3.066479806614716e+01
        a6 :: a
a6 =  2.506628277459239e+00

        b1 :: a
b1 = -5.447609879822406e+01
        b2 :: a
b2 =  1.615858368580409e+02
        b3 :: a
b3 = -1.556989798598866e+02
        b4 :: a
b4 =  6.680131188771972e+01
        b5 :: a
b5 = -1.328068155288572e+01

        c1 :: a
c1 = -7.784894002430293e-03
        c2 :: a
c2 = -3.223964580411365e-01
        c3 :: a
c3 = -2.400758277161838e+00
        c4 :: a
c4 = -2.549732539343734e+00
        c5 :: a
c5 =  4.374664141464968e+00
        c6 :: a
c6 =  2.938163982698783e+00

        d1 :: a
d1 =  7.784695709041462e-03
        d2 :: a
d2 =  3.224671290700398e-01
        d3 :: a
d3 =  2.445134137142996e+00
        d4 :: a
d4 =  3.754408661907416e+00

        pLow :: a
pLow = 0.02425

        nan :: a
nan = 0a -> a -> a
forall a. Fractional a => a -> a -> a
/0

    in  if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0 then
            a
nan
        else if a
p a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then
            -1a -> a -> a
forall a. Fractional a => a -> a -> a
/0
        else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
pLow then
            let q :: a
q = a -> a
forall a. Floating a => a -> a
sqrt(-2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
log a
p)
            in  (((((a
c1a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c2)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c3)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c4)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c5)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
c6) a -> a -> a
forall a. Fractional a => a -> a -> a
/
                 ((((a
d1a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
d2)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
d3)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+a
d4)a -> a -> a
forall a. Num a => a -> a -> a
*a
qa -> a -> a
forall a. Num a => a -> a -> a
+1)
        else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 1 a -> a -> a
forall a. Num a => a -> a -> a
- a
pLow then
            let q :: a
q = a
p a -> a -> a
forall a. Num a => a -> a -> a
- 0.5
                r :: a
r = a
qa -> a -> a
forall a. Num a => a -> a -> a
*a
q
            in  (((((a
a1a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a2)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a3)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a4)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a5)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
a6)a -> a -> a
forall a. Num a => a -> a -> a
*a
q a -> a -> a
forall a. Fractional a => a -> a -> a
/
                (((((a
b1a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b2)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b3)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b4)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+a
b5)a -> a -> a
forall a. Num a => a -> a -> a
*a
ra -> a -> a
forall a. Num a => a -> a -> a
+1)
        else if a
p a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 then
            - a -> a
forall a. (Ord a, Floating a) => a -> a
invnormcdf (1 a -> a -> a
forall a. Num a => a -> a -> a
- a
p)
        else
            a
nan


-- | Try to estimate complexity of a whole from a sample.  Suppose we
-- sampled @total@ things and among those @singles@ occured only once.
-- How many different things are there?
--
-- Let the total number be @m@.  The copy number follows a Poisson
-- distribution with paramter @\lambda@.  Let \( z := e^{\lambda} \), then
-- we have:
--
-- \[
--   P( 0 ) = e^{-\lambda} = \frac{1}{z}                    \\
--   P( 1 ) = \lambda e^{-\lambda} = \frac{\ln z}{z}                   \\
--   P(\ge 1) = 1 - e^{-\lambda} = 1 - \frac{1}{z}                    \\
-- \]
-- \[
--   \mbox{singles} = m \frac{\ln z}{z}                   \\
--   \mbox{total}   = m \left( 1 - \frac{1}{z} \right)                  \\
-- \]
-- \[
--   D := \frac{\mbox{total}}{\mbox{singles}} = (1 - \frac{1}{z}) * \frac{z}{\ln z}                  \\
--   f := z - 1 - D \ln z = 0
-- \]
--
-- To get @z@, we solve using Newton iteration and then substitute to
-- get @m@:
--
-- \[
--   df/dz = 1 - D/z                                    \\
--   z' = z - \frac{ z (z - 1 - D \ln z) }{ z - D }     \\
--   m = \mbox{singles} * \frac{z}{\ln z}
-- \]
--
-- It converges as long as the initial @z@ is large enough, and @10D@
-- (in the line for @zz@ below) appears to work well.

estimateComplexity :: (Integral a, Floating b, Ord b) => a -> a -> Maybe b
estimateComplexity :: a -> a -> Maybe b
estimateComplexity total :: a
total singles :: a
singles | a
total   a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
singles = Maybe b
forall a. Maybe a
Nothing
                                 | a
singles a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= 0       = Maybe b
forall a. Maybe a
Nothing
                                 | Bool
otherwise          = b -> Maybe b
forall a. a -> Maybe a
Just b
m
  where
    d :: b
d = a -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
total b -> b -> b
forall a. Fractional a => a -> a -> a
/ a -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
singles
    step :: b -> b
step z :: b
z = b
z b -> b -> b
forall a. Num a => a -> a -> a
* (b
z b -> b -> b
forall a. Num a => a -> a -> a
- 1 b -> b -> b
forall a. Num a => a -> a -> a
- b
d b -> b -> b
forall a. Num a => a -> a -> a
* b -> b
forall a. Floating a => a -> a
log b
z) b -> b -> b
forall a. Fractional a => a -> a -> a
/ (b
z b -> b -> b
forall a. Num a => a -> a -> a
- b
d)
    iter :: b -> b
iter z :: b
z = case b -> b
step b
z of zd :: b
zd | b -> b
forall a. Num a => a -> a
abs b
zd b -> b -> Bool
forall a. Ord a => a -> a -> Bool
< 1e-12 -> b
z
                               | Bool
otherwise -> b -> b
iter (b -> b) -> b -> b
forall a b. (a -> b) -> a -> b
$! b
zb -> b -> b
forall a. Num a => a -> a -> a
-b
zd
    zz :: b
zz = b -> b
iter (b -> b) -> b -> b
forall a b. (a -> b) -> a -> b
$! 10b -> b -> b
forall a. Num a => a -> a -> a
*b
d
    m :: b
m = a -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
singles b -> b -> b
forall a. Num a => a -> a -> a
* b
zz b -> b -> b
forall a. Fractional a => a -> a -> a
/ b -> b
forall a. Floating a => a -> a
log b
zz


-- | Computes \( \ln \left( e^x + e^y \right) \) without leaving the log domain and
-- hence without losing precision.
infixl 5 <#>
{-# INLINE (<#>) #-}
(<#>) :: (Floating a, Ord a) => a -> a -> a
x :: a
x <#> :: a -> a -> a
<#> y :: a
y = if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y then a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. (Floating a, Ord a) => a -> a
log1pexp (a
ya -> a -> a
forall a. Num a => a -> a -> a
-a
x) else a
y a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. (Floating a, Ord a) => a -> a
log1pexp (a
xa -> a -> a
forall a. Num a => a -> a -> a
-a
y)

-- | Computes @log (1+x)@ to a relative precision of @10^-8@ even for
-- very small @x@.  Stolen from <http://www.johndcook.com/cpp_log_one_plus_x.html>
{-# INLINE log1p #-}
log1p :: (Floating a, Ord a) => a -> a
log1p :: a -> a
log1p x :: a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< -1 = String -> a
forall a. HasCallStack => String -> a
error "log1p: argument must be greater than -1"
        -- x is large enough that the obvious evaluation is OK:
        | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> 0.0001 Bool -> Bool -> Bool
|| a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< -0.0001 = a -> a
forall a. Floating a => a -> a
log (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ 1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
x
        -- Use Taylor approx. log(1 + x) = x - x^2/2 with error roughly x^3/3
        -- Since |x| < 10^-4, |x|^3 < 10^-12, relative error less than 10^-8:
        | Bool
otherwise = (1 a -> a -> a
forall a. Num a => a -> a -> a
- 0.5a -> a -> a
forall a. Num a => a -> a -> a
*a
x) a -> a -> a
forall a. Num a => a -> a -> a
* a
x


-- | Computes \( e^x - 1 \) to a relative precision of @10^-10@ even for
-- very small @x@.  Stolen from <http://www.johndcook.com/cpp_expm1.html>
{-# INLINE expm1 #-}
expm1 :: (Floating a, Ord a) => a -> a
expm1 :: a -> a
expm1 x :: a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> -0.00001 Bool -> Bool -> Bool
&& a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0.00001 = (1 a -> a -> a
forall a. Num a => a -> a -> a
+ 0.5 a -> a -> a
forall a. Num a => a -> a -> a
* a
x) a -> a -> a
forall a. Num a => a -> a -> a
* a
x       -- Taylor approx
        | Bool
otherwise                   = a -> a
forall a. Floating a => a -> a
exp a
x a -> a -> a
forall a. Num a => a -> a -> a
- 1               -- direct eval

-- | Computes \( \ln (1 - e^x) \), following Martin Mächler.
{-# INLINE log1mexp #-}
log1mexp :: (Floating a, Ord a) => a -> a
log1mexp :: a -> a
log1mexp x :: a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> - a -> a
forall a. Floating a => a -> a
log 2 = a -> a
forall a. Floating a => a -> a
log (- a -> a
forall a. (Floating a, Ord a) => a -> a
expm1 a
x)
           | Bool
otherwise   = a -> a
forall a. (Floating a, Ord a) => a -> a
log1p (- a -> a
forall a. Floating a => a -> a
exp a
x)

-- | Computes \( \ln (1 + e^x) \), following Martin Mächler.
{-# INLINE log1pexp #-}
log1pexp :: (Floating a, Ord a) => a -> a
log1pexp :: a -> a
log1pexp x :: a
x | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<=  -37 = a -> a
forall a. Floating a => a -> a
exp a
x
           | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<=   18 = a -> a
forall a. (Floating a, Ord a) => a -> a
log1p (a -> a) -> a -> a
forall a b. (a -> b) -> a -> b
$ a -> a
forall a. Floating a => a -> a
exp a
x
           | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= 33.3 = a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
exp (-a
x)
           | Bool
otherwise = a
x


-- | Computes \( \ln ( \sum_i e^{x_i} ) \) sensibly.  The list must be
-- sorted in descending(!) order.
{-# INLINE lsum #-}
lsum :: (Floating a, Ord a) => [a] -> a
lsum :: [a] -> a
lsum = (a -> a -> a) -> [a] -> a
forall a. (a -> a -> a) -> [a] -> a
foldl1' (\x :: a
x y :: a
y -> if a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y then a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. (Floating a, Ord a) => a -> a
log1pexp (a
ya -> a -> a
forall a. Num a => a -> a -> a
-a
x) else a
forall a. a
err)
    where err :: a
err = String -> a
forall a. HasCallStack => String -> a
error "lsum: argument list must be in descending order"

-- | Computes \( \ln \left( c e^x + (1-c) e^y \right) \).
{-# INLINE llerp #-}
llerp :: (Floating a, Ord a) => a -> a -> a -> a
llerp :: a -> a -> a -> a
llerp c :: a
c x :: a
x y :: a
y | a
c a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= 0.0  = a
y
            | a
c a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= 1.0  = a
x
            | a
x a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
y    = a -> a
forall a. Floating a => a -> a
log     a
c  a -> a -> a
forall a. Num a => a -> a -> a
+ a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. (Floating a, Ord a) => a -> a
log1p ( (1a -> a -> a
forall a. Num a => a -> a -> a
-a
c)a -> a -> a
forall a. Fractional a => a -> a -> a
/a
c a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
exp (a
ya -> a -> a
forall a. Num a => a -> a -> a
-a
x) )        -- Hmm.
            | Bool
otherwise = a -> a
forall a. (Floating a, Ord a) => a -> a
log1p (-a
c) a -> a -> a
forall a. Num a => a -> a -> a
+ a
y a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. (Floating a, Ord a) => a -> a
log1p ( a
ca -> a -> a
forall a. Fractional a => a -> a -> a
/(1a -> a -> a
forall a. Num a => a -> a -> a
-a
c) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
exp (a
xa -> a -> a
forall a. Num a => a -> a -> a
-a
y) )        -- Hmm.

-- | Binomial coefficient: \( \mbox{choose n k} = \frac{n!}{(n-k)! k!} \)
{-# INLINE choose #-}
choose :: Integral a => a -> a -> a
choose :: a -> a -> a
choose n :: a
n k :: a
k = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [a
na -> a -> a
forall a. Num a => a -> a -> a
-a
ka -> a -> a
forall a. Num a => a -> a -> a
+1 .. a
n] a -> a -> a
forall a. Integral a => a -> a -> a
`div` [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [2..a
k]