diversity-0.8.1.0: src/src-lib/Math/Diversity/Diversity.hs
-- Diversity module.
-- By G.W. Schwartz
--
{- | Collection of functions pertaining to finding the diversity of samples.
-}
{-# LANGUAGE BangPatterns #-}
module Math.Diversity.Diversity ( hamming
, richness
, diversity
, diversityOfMap
, diversityOfMapWithContribution
, chao1
, chao2
, chao1Var
, chao2Var
, rarefactionCurve
, rarefactionSampleCurve
, rarefactionViable
, individualG
, sampleG
, minRarefaction ) where
-- Built-in
import Data.List
import qualified Data.Set as Set
import qualified Data.Map.Strict as Map
import Numeric.SpecFunctions (choose)
-- Cabal
import qualified Data.List.Ordered as LO
-- Local
import Math.Diversity.RandomSampling
import Math.Diversity.Types
-- | Takes two strings, returns Hamming distance
hamming :: String -> String -> Int
hamming xs ys = length $ filter not $ zipWith (==) xs ys
-- | Fast product division
productDivision :: Double -> [Integer] -> [Integer] -> Double
productDivision acc [] [] = acc
productDivision acc [] (y:ys) = (acc / fromInteger y)
* productDivision acc [] ys
productDivision acc (x:xs) [] = acc * fromInteger x * productDivision acc xs []
productDivision acc (x:xs) (y:ys)
| x == y = productDivision acc xs ys
| otherwise = (fromInteger x / fromInteger y) * productDivision acc xs ys
-- | Returns the diversity of a list of things
diversity :: (Ord a) => Double -> [a] -> Double
diversity order sample
| null sample = 0
| order == 1 = exp . h . speciesList $ sample
| otherwise = ( Map.foldl' (+) 0
. Map.map ((** order) . p_i)
. speciesList
$ sample )
** pow
where
pow = 1 / (1 - order)
h = negate
. Map.foldl' (+) 0
. Map.map (\x -> p_i x * log (p_i x))
p_i x = (x :: Double) /
((Map.foldl' (+) 0 . speciesList $ sample) :: Double)
speciesList = Map.fromListWith (+) . map (\x -> (x, 1))
-- | Returns the diversity of a map of the species and how many times it
-- appears
diversityOfMap :: (Ord a) => Double -> Map.Map a Int -> Double
diversityOfMap order sample
| Map.null sample = 0
| order == 1 = exp . h $ sample
| otherwise = ( Map.foldl' (+) 0
. Map.map ((** order) . p_i)
$ sample )
** pow
where
pow = 1 / (1 - order)
h = negate
. Map.foldl' (+) 0
. Map.map ( \x -> p_i x * log (p_i x))
p_i x = fromIntegral x / fromIntegral (Map.foldl' (+) 0 sample)
-- | Returns the diversity of a map of the species and how many times it
-- appears, along with the contribution of each species.
diversityOfMapWithContribution
:: (Ord a)
=> Double -> Map.Map a Int -> (Double, Map.Map a Double)
diversityOfMapWithContribution order sample
| Map.null sample = (0, Map.empty)
| order == 1 = (exp . h $ sample, contrib)
| otherwise = ( ( Map.foldl' (+) 0
. Map.map ((** order) . p_i)
$ sample
)
** pow
, contrib
)
where
contrib = Map.map
( (/ (fromIntegral $ Map.foldl' (+) 0 sample))
. fromIntegral
)
sample
pow = 1 / (1 - order)
h = negate
. Map.foldl' (+) 0
. Map.map ( \x -> p_i x * log (p_i x))
p_i x = fromIntegral x / fromIntegral (Map.foldl' (+) 0 sample)
-- | Returns the richness of the observed data
richness :: (Ord a, Ord b) => Map.Map (a, b) c -> Int
richness = Map.size . Map.mapKeys snd
-- | Returns the map of species with how many samples they appear in
overlapSampleMap :: (Ord a, Ord b) => Map.Map (a, b) Int -> Map.Map b Int
overlapSampleMap = Map.mapKeysWith (+) snd . Map.map (const 1)
-- | Returns the number of a that appear x times
abundanceFreq :: (Ord a) => Int -> Map.Map a Int -> Int
abundanceFreq x = Map.size . Map.filter (== x)
-- | Returns the number of b that appear in x number of a. Notice that this
-- function takes in a normal frequency map, as it converts it with
-- overlapSampleMap.
overlapFreq :: (Ord a, Ord b) => Int -> Map.Map (a, b) Int -> Int
overlapFreq x = Map.size . Map.filter (== x) . overlapSampleMap
-- | Returns the chao1 estimator of a map of the species and how many times
-- it appears
chao1 :: (Ord a) => Map.Map a Int -> Double
chao1 sample
| f2 > 0 = (f1 ** 2) / (2 * f2)
| otherwise = (f1 * (f1 - 1)) / (2 * (f2 + 1))
where
f1 = fromIntegral . abundanceFreq 1 $ sample
f2 = fromIntegral . abundanceFreq 2 $ sample
-- | Returns the chao2 estimator of a map of the sample labeled species
-- (sample, species) and how many times it appears. This will calculate the
-- overlap for you, so if you don't have the number of times it appears it
-- does not matter, you can set it to 1 and get the same result as it's all
-- about overlapping samples.
chao2 :: (Ord a, Ord b) => Map.Map (a, b) Int -> Double
chao2 samples
| q2 > 0 = ((t - 1) / t) * ((q1 ** 2) / (2 * q2))
| otherwise = ((t - 1) / t) * ((q1 * (q1 - 1)) / (2 * (q2 + 1)))
where
q1 = fromIntegral . overlapFreq 1 $ samples
q2 = fromIntegral . overlapFreq 2 $ samples
-- Saves time so don't have to recalculate
t = fromIntegral . Map.size . Map.mapKeys fst $ samples
-- | Returns the chao1 estimator variance of a map of the species
-- and how many times each one appears.
chao1Var :: (Ord a) => Map.Map a Int -> Double
chao1Var sample
| f2 > 0 = f2
* ( ((1 / 2) * ((n - 1) / n) * ((f1 / f2) ** 2))
+ ((((n - 1) / n) ** 2) * ((f1 / f2) ** 3))
+ ((1 / 4) * (((n - 1) / n) ** 2) * ((f1 / f2) ** 4))
)
| otherwise = (((n - 1) / n) * ((f1 * (f1 - 1)) / 2))
+ ((((n - 1) / n) ** 2) * ((f1 * (((2 * f1) - 1) ** 2)) / 4))
+ ((((n - 1) / n) ** 2) * ((f1 ** 4) / (4 * sest)))
where
sest = fromIntegral (Map.size sample) + chao1 sample
f1 = fromIntegral . abundanceFreq 1 $ sample
f2 = fromIntegral . abundanceFreq 2 $ sample
-- Saves time so don't have to recalculate
n = fromIntegral . Map.foldl' (+) 0 $ sample
-- | Returns the chao2 estimator variance of a map of the sample labeled species
-- (sample, species) and how many times it appears.
chao2Var :: (Ord a, Ord b) => Map.Map (a, b) Int -> Double
chao2Var samples
| q2 > 0 = q2
* ( ((1 / 2) * ((t - 1) / t) * ((q1 / q2) ** 2))
+ ((((t - 1) / t) ** 2) * ((q1 / q2) ** 3))
+ ((1 / 4) * ((t - 1 / t) ** 2) * ((q1 / q2) ** 4))
)
| otherwise = (((t - 1) / t) * ((q1 * (q1 - 1)) / 2))
+ ((((t - 1) / t) ** 2) * ((q1 * (((2 * q1) - 1) ** 2)) / 4))
+ ((((t - 1) / t) ** 2) * ((q1 ** 4) / (4 * sest)))
where
sest = fromIntegral (richness samples) + chao2 samples
q1 = fromIntegral . overlapFreq 1 $ samples
q2 = fromIntegral . overlapFreq 2 $ samples
-- Saves time so don't have to recalculate
t = fromIntegral . Map.size . Map.mapKeys fst $ samples
-- | Binomial for small or large numbers (slow but works for big numbers,
-- fast but works for small numbers)
specialBinomial :: Bool -> Integer -> Integer -> Integer -> Double
specialBinomial False n_total g n = productDivision 1 num den
where
num = [(n_total - g - n + 1)..(n_total - g)]
den = [(n_total - n + 1)..n_total]
specialBinomial True n_total g n = choose
(fromIntegral n_total - fromIntegral g)
(fromIntegral n)
-- | Returns the rarefaction curve for each position in a list
rarefactionCurve :: Bool
-> Int
-> Integer
-> Integer
-> Integer
-> Map.Map (Sample, Fragment) Int
-> IO [(Int, Maybe (Double, Double))]
rarefactionCurve !fastBin !runs !start !interval !end !sample =
mapM rarefact
. LO.nubSort
$ n_total : [start,(start + interval)..finish]
where
rarefact !n
| n == 0 = return (fromIntegral n, Just (0, 0))
| n == 1 = return (fromIntegral n, Just (1, 0))
| n == n_total = return (fromIntegral n, Just (k, 0))
| n > n_total = return (fromIntegral n, Just (estimation n, 0))
| runs == 0 = return (fromIntegral n, Just (k - inner n, 0))
| otherwise = do -- Empirical version
statTuple <- subsampleES
runs
(fromIntegral n_total)
(fromIntegral n)
. concatMap snd
. Map.toAscList
. Map.mapWithKey (\(_, f) x -> replicate x f)
$ sample
return (fromIntegral n, statTuple)
inner n = ( \x -> if fastBin
then x / choose (fromIntegral n_total) (fromIntegral n)
else x )
. sum
. map (\g -> specialBinomial fastBin n_total (fromIntegral g) n)
$ grouped
-- Unreadable unless I break the 80 column rule
estimation n = fromIntegral n_total + (chao1 sample * (1 - exp (((fromIntegral n - fromIntegral n_total) / (- fromIntegral n_total)) * (fromIntegral (abundanceFreq 1 sample) / chao1 sample))))
finish = if end == 0 then n_total else end
n_total = fromIntegral . Map.foldl' (+) 0 $ sample
k = fromIntegral . Map.size $ sample
grouped = Map.elems sample
-- | Each sample has a collection of species, return a list of these maps
getSampleContents :: Map.Map (Sample, Fragment) Int -> [Set.Set Fragment]
getSampleContents = Map.elems
. Map.fromListWith Set.union
. map (\(x, y) -> (x, Set.singleton y))
. Map.keys
-- | Returns the rarefaction curve for each position in a list
rarefactionSampleCurve :: Bool
-> Int
-> Int
-> Int
-> Map.Map (Sample, Fragment) Int
-> IO [(Int, Maybe (Double, Double))]
rarefactionSampleCurve !fastBin !start !interval !end !ls =
mapM rarefact
. LO.nubSort
$ t_total : [start,(start + interval)..finish]
where
rarefact !t
| t == 0 = return (t, Just (0, 0))
| t == t_total = return (t, Just (sobs, 0))
| t > t_total = return (t, Just (estimation t, 0))
| otherwise = return (t, Just (sobs - inner t, 0))
inner t = ( \x -> if fastBin
then x / choose t_total t
else x )
. sum
. map ( \s -> specialBinomial
fastBin
(fromIntegral t_total)
(numHave s samples)
(fromIntegral t) )
$ speciesList
-- Unreadable unless I break the 80 column rule
estimation t = sobs + (chao2 ls * (1 - exp ((- (fromIntegral t - fromIntegral t_total) * fromIntegral (overlapFreq 1 ls)) / (fromIntegral (overlapFreq 1 ls) + (fromIntegral t_total * chao2 ls)))))
finish = if end == 0 then t_total else end
numHave s = genericLength . filter (Set.member s)
sobs = fromIntegral $ richness ls
speciesList = Map.keys . Map.mapKeys snd $ ls
t_total = genericLength samples
samples = getSampleContents ls
-- | Calculates the percent of the curve that is above 95% of height of the
-- curve
rarefactionViable :: [Double] -> Double
rarefactionViable !xs = (genericLength valid / genericLength xs) * 100
where
valid = dropWhile (< (0.95 * last xs)) xs
-- | Returns the number of individuals needed to get the proportion g of
-- the estimated total richness of the assemblage. Sobs / Sest < g < 1
individualG :: Double -> Map.Map (Sample, Fragment) Int -> Double
individualG g sample
| sobs / sest >= g = 0
| otherwise = ((sobs * f1) / (2 * f2)) * log (f0 / ((1 - g) * sest))
where
sest = sobs + f0
sobs = fromIntegral $ richness sample
f2 = fromIntegral . abundanceFreq 2 $ sample
f1 = fromIntegral . abundanceFreq 1 $ sample
f0 = chao1 sample
-- | Returns the number of samples needed to get the proportion g of
-- the estimated total richness of the assemblage. Sobs / Sest < g < 1
sampleG :: Double -> Map.Map (Sample, Fragment) Int -> Double
sampleG g sample
| sobs / sest >= g = 0
| otherwise = log (1 - ((t / (t - 1)) * ((2 * q2) / (q1 ** 2)) * ((g * sest) - sobs)))
/ log (1 - ((2 * q2) / (((t - 1) * q1) + (2 * q2))))
where
t = fromIntegral . Map.size . Map.mapKeys fst $ sample
sest = sobs + q0
sobs = fromIntegral $ richness sample
q2 = fromIntegral . overlapFreq 2 $ sample
q1 = fromIntegral . overlapFreq 1 $ sample
q0 = chao2 sample
-- | Returns the number of samples needed before a new sample returns less
-- than x new species. Warning, goes forever until threshold is met!
minRarefaction :: Bool
-> Bool
-> Int
-> Map.Map (Sample, Fragment) Int
-> Double
-> Int
-> IO Int
minRarefaction _ _ (-1) _ _ _ = return (-1)
minRarefaction bySample fastBin threshold sample !oldRare !count = do
newRare <- fmap (maybe (-1) fst . snd . head)
. rarefaction bySample count
$ sample
if newRare == (-1)
then return (-1)
else if newRare - oldRare < fromIntegral threshold
then return count
else minRarefaction True fastBin threshold sample newRare (count + 1)
where
rarefaction True x = rarefactionSampleCurve fastBin x 1 x
rarefaction False x = rarefactionCurve
fastBin
0
(fromIntegral x)
1
(fromIntegral x)