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 )
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
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
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
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)
{-# 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"
| 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
| 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
{-# 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
| Bool
otherwise = a -> a
forall a. Floating a => a -> a
exp a
x a -> a -> a
forall a. Num a => a -> a -> a
- 1
{-# 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)
{-# 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
{-# 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"
{-# 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) )
| 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) )
{-# 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]