module Pareto
( controlLimit
, entropy
, ratio
, pareto
, causeEffect
, causesMaxConcentration
, effectsMaxConcentration
) where

import Data.List (sort)

-- | Determine control limit.
--
-- >>> controlLimit 10
-- 2.7709505944546686
controlLimit x = -1 * majority - minority
  where
    factor   = x / 5
    majority = 0.6 * (logBase 2 (0.6 / factor))
    minority = 0.4 * (logBase 2 (0.4 / (4 * factor)))

-- | Determine entropy for list of numbers
--
-- >>> entropy [0.3, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05] 
-- 2.770950594454668
--
-- >>> entropy [789, 621, 109, 65, 45, 30, 27, 15, 12, 9]
-- 1.9593816735406657
entropy numbers = foldl (\acc x -> acc - subtrahend x) 0 numbers
  where
    s = sum numbers
    subtrahend x = (x/s) * logBase 2 (x/s)

-- | entropy divided by controlLimit.
--
-- >>> ratio [0.3, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05]
-- 1.0
ratio numbers = etr / cl
  where
    etr = roundTo 6 $ entropy numbers
    cl  = roundTo 6 $ controlLimit $ fromIntegral $ length numbers 

-- | Returns True if pareto distribution is present otherwise False.
--
-- >>> pareto [0.3, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05]
-- True
pareto numbers = ratio numbers <= 1

-- | Determine effects for specified ratio of causes.
--
-- >>> effects 0.2 [0.3, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05]
-- 0.6
effects causes numbers = last $ take halt eff
  where
    caus = makeCauses $ fromIntegral $ length numbers
    eff  = makeEffects numbers
    halt = l - (fromIntegral $ length $ dropWhile (<= causes) caus)
    l :: Num a => a
    l    = fromIntegral $ length numbers

-- | Determine causes for specified ratio of effects.
--
-- >>> causes 0.6 [0.3, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05]
-- 0.2
causes effects numbers = last $ take halt caus
  where
    caus = makeCauses $ fromIntegral $ length numbers
    eff  = makeEffects numbers
    halt = fromIntegral $ length $ takeWhile (<= effects) eff
    l :: Num a => a
    l    = fromIntegral $ length numbers

-- | Return a list of causes of specified length.
--
-- >>> makeCauses 10
-- [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
makeCauses :: (Fractional b, Enum b) => b -> [b]
makeCauses l = map (/l) [1..l]

-- | Return a list of sorted cumsum in proportion to sum of list of numbers.
--
-- >>> makeEffects  [0.3, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05]
-- [0.3,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0]
makeEffects numbers = map (roundTo 3 . (/s)) cumsum
  where
    s      = sum numbers
    sorted = reverse $ sort numbers
    cumsum = scanl1 (+) sorted

-- | Which causes have the maximum concentration (rank * value)?
--
-- >>> causesMaxConcentration [789, 621, 109, 65, 45, 30, 27, 15, 12, 9]
-- 0.2
causesMaxConcentration :: [Float] -> Float
causesMaxConcentration numbers = ordinal / (fromIntegral $ length numbers)
  where
    concentration = makeConcentration numbers
    peak          = maximum concentration
    equalMax      = map (peak ==) concentration
    ordinal       = fromIntegral $ (length $ takeWhile (== False) equalMax) + 1

-- | Which effects have the maximum concentration (rank * value)?
--
-- >>> effectsMaxConcentration [789, 621, 109, 65, 45, 30, 27, 15, 12, 9]
-- 0.819
effectsMaxConcentration :: Num a => [Float] -> Float
effectsMaxConcentration numbers = effects (causesMaxConcentration numbers) numbers

-- | Return list of cause-effect pairs.
--
-- >>> causeEffect [789, 621, 109, 65, 45, 30, 27, 15, 12, 9]
-- [(0.1,0.458),(0.2,0.819),(0.3,0.882),(0.4,0.92),(0.5,0.946),(0.6,0.963),(0.7,0.979),(0.8,0.988),(0.9,0.995),(1.0,1.0)]
causeEffect :: (RealFrac a1, Fractional b, Fractional a, Enum a) => [a1] -> [(a, b)]
causeEffect numbers = zip caus eff
  where
    caus = makeCauses l
    eff  = makeEffects numbers
    l    = fromIntegral $ length numbers

-- | Calculate concentration (rank * value) for list of numbers.
--
-- >>> makeConcentration [789, 621, 109, 65, 45, 30, 27, 15, 12, 9]
-- [789,1242,327,260,225,180,189,120,108,90]
makeConcentration numbers = zipWith (*) [1..l] sorted
  where
    sorted = reverse $ sort numbers
    l = fromIntegral $ length numbers

-- Helpers
    
-- | Round number to specified number of decimals.
--
-- >>> roundTo 6 2.7709505944546686 
-- 2.770951
roundTo decimals number = fromIntegral (round $ number * 10^decimals) / 10^decimals