module Bio.Util.Numeric (
wilson, invnormcdf, choose,
estimateComplexity, showNum, showOOM,
log1p, expm1, (<#>),
lsum, llerp,
sigmoid2, isigmoid2
) where
import Data.List ( foldl1' )
import Data.Char ( intToDigit )
wilson :: Double -> Int -> Int -> (Double, Double, Double)
wilson c x n = ( (m h) / d, p, (m + h) / d )
where
nn = fromIntegral n
p = fromIntegral x / nn
z = invnormcdf (1c*0.5)
h = z * sqrt (( p * (1p) + 0.25*z*z / nn ) / nn)
m = p + 0.5 * z * z / nn
d = 1 + z * z / nn
showNum :: Show a => a -> String
showNum = triplets [] . reverse . show
where
triplets acc [] = acc
triplets acc (a:[]) = a:acc
triplets acc (a:b:[]) = b:a:acc
triplets acc (a:b:c:[]) = c:b:a:acc
triplets acc (a:b:c:s) = triplets (',':c:b:a:acc) s
showOOM :: Double -> String
showOOM x | x < 0 = '-' : showOOM (negate x)
| otherwise = findSuffix (x*10) ".kMGTPEZY"
where
findSuffix _ [] = "many"
findSuffix y (s:ss) | y < 100 = intToDigit (round y `div` 10) : case (round y `mod` 10, s) of
(0,'.') -> [] ; (0,_) -> [s] ; (d,_) -> [s, intToDigit d]
| y < 1000 = intToDigit (round y `div` 100) : intToDigit ((round y `mod` 100) `div` 10) :
if s == '.' then [] else [s]
| y < 10000 = intToDigit (round y `div` 1000) : intToDigit ((round y `mod` 1000) `div` 100) :
'0' : if s == '.' then [] else [s]
| otherwise = findSuffix (y*0.001) ss
invnormcdf :: (Ord a, Floating a) => a -> a
invnormcdf p =
let a1 = 3.969683028665376e+01
a2 = 2.209460984245205e+02
a3 = 2.759285104469687e+02
a4 = 1.383577518672690e+02
a5 = 3.066479806614716e+01
a6 = 2.506628277459239e+00
b1 = 5.447609879822406e+01
b2 = 1.615858368580409e+02
b3 = 1.556989798598866e+02
b4 = 6.680131188771972e+01
b5 = 1.328068155288572e+01
c1 = 7.784894002430293e-03
c2 = 3.223964580411365e-01
c3 = 2.400758277161838e+00
c4 = 2.549732539343734e+00
c5 = 4.374664141464968e+00
c6 = 2.938163982698783e+00
d1 = 7.784695709041462e-03
d2 = 3.224671290700398e-01
d3 = 2.445134137142996e+00
d4 = 3.754408661907416e+00
pLow = 0.02425
nan = 0/0
in if p < 0 then
nan
else if p == 0 then
1/0
else if p < pLow then
let q = sqrt(2 * log p)
in (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) /
((((d1*q+d2)*q+d3)*q+d4)*q+1)
else if p < 1 pLow then
let q = p 0.5
r = q*q
in (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q /
(((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1)
else if p <= 1 then
invnormcdf (1 p)
else
nan
estimateComplexity :: (Integral a, Floating b, Ord b) => a -> a -> Maybe b
estimateComplexity total singles | total <= singles = Nothing
| singles <= 0 = Nothing
| otherwise = Just m
where
d = fromIntegral total / fromIntegral singles
step z = z * (z 1 d * log z) / (z d)
iter z = case step z of zd | abs zd < 1e-12 -> z
| otherwise -> iter $! zzd
zz = iter $! 10*d
m = fromIntegral singles * zz / log zz
infixl 5 <#>
(<#>) :: (Floating a, Ord a) => a -> a -> a
x <#> y = if x >= y then x + log1p (exp (yx)) else y + log1p (exp (xy))
log1p :: (Floating a, Ord a) => a -> a
log1p x | x < 1 = error "log1p: argument must be greater than -1"
| x > 0.0001 || x < 0.0001 = log $ 1 + x
| otherwise = (1 0.5*x) * x
expm1 :: (Floating a, Ord a) => a -> a
expm1 x | x > 0.00001 && x < 0.00001 = (1 + 0.5 * x) * x
| otherwise = exp x 1
lsum :: (Floating a, Ord a) => [a] -> a
lsum xs = foldl1' (\x y -> if x >= y then x + log1p (exp (yx)) else err) xs
where err = error $ "lsum: argument list must be in descending order"
llerp :: (Floating a, Ord a) => a -> a -> a -> a
llerp c x y | c <= 0.0 = y
| c >= 1.0 = x
| x >= y = log c + x + log1p ( (1c)/c * exp (yx) )
| otherwise = log1p (c) + y + log1p ( c/(1c) * exp (xy) )
choose :: Integral a => a -> a -> a
n `choose` k = product [nk+1 .. n] `div` product [2..k]
sigmoid2 :: (Num a, Fractional a, Floating a) => a -> a
sigmoid2 x = y*y where y = (exp x 1) / (exp x + 1)
isigmoid2 :: (Num a, Fractional a, Floating a) => a -> a
isigmoid2 y = log $ (1 + sqrt y) / (1 sqrt y)