-- |
-- Module     : Simulation.Aivika.Experiment.Histogram
-- Copyright  : Copyright (c) 2012-2017, David Sorokin <david.sorokin@gmail.com>
-- License    : BSD3
-- Maintainer : David Sorokin <david.sorokin@gmail.com>
-- Stability  : experimental
-- Tested with: GHC 8.0.1
--
-- This module computes the histogram by the 
-- specified data and strategy applied for such computing.
--
-- The code in this module is essentially based on the 
-- <http://hackage.haskell.org/package/Histogram> package
-- by Mike Izbicki, who kindly agreed to re-license 
-- his library under BSD3, which allowed me to use his 
-- code and comments with some modifications. 

module Simulation.Aivika.Experiment.Histogram 
    (-- * Creating Histograms
     Histogram(..), 
     histogram, 
     histogramBinSize, 
     histogramNumBins,
     -- * Binning Strategies
     BinningStrategy(..),
     binSturges, 
     binDoane, 
     binSqrt, 
     binScott) where

import Data.List
import Data.Monoid
import qualified Data.Map as Map

import Numeric

-------------------------------------------------------------------------------

-- | Holds all the information needed to plot the histogram 
-- for a list of different series. Each series produces its 
-- own item in the resuling @[Int]@ list that may contain
-- zeros.
type Histogram = [(Double, [Int])]

-------------------------------------------------------------------------------
-- Bin counters; check out http://en.wikipedia.org/wiki/Histogram#Number_of_bins_and_width

-- | The strategy applied to calculate the histogram bins.
type BinningStrategy = [Double] -> Int

-- | This is only a helper function to convert strategies that 
-- specify bin width into strategies that specify the number of bins.
stratFromBinWidth :: [Double] -> Double -> Int
stratFromBinWidth xs h = ceiling $ (maximum xs - minimum xs) / h

-- | Sturges' binning strategy is the least computational work, 
-- but recommended for only normal data.
binSturges :: BinningStrategy
binSturges xs = ceiling $ logBase 2 n + 1
    where n = fromIntegral $ length xs
          
-- | Doane's binning strategy extends Sturges' for non-normal data.  
-- It takes a little more time because it must calculate the kurtosis 
-- (peakkiness) of the distribution.
binDoane :: BinningStrategy
binDoane xs = ceiling $ 1 + log n + log (1 + a * ((n/6)**(1/2)))
    where 
        n = fromIntegral $ length xs
        a = kurtosis xs
        
-- | Using the sqrt of the number of samples is not supported by any 
-- theory, but is commonly used by excel and other histogram making software.
binSqrt :: BinningStrategy
binSqrt xs = round $ sqrt n
    where 
        n = fromIntegral $ length xs

-- | Scott's rule is the optimal solution for normal data, but requires 
-- more computation than Sturges'.
binScott :: BinningStrategy
binScott xs = stratFromBinWidth xs $ 3.5 * stddev xs / (n**(1/3))
    where 
        n = fromIntegral $ length xs
        
-------------------------------------------------------------------------------
-- create the histogram data

-- | Creates a histogram by specifying the list of series.  Call it with one of 
-- the binning strategies that is appropriate to the type of data you have.  
-- If you don't know, then try using 'binSturges'.
histogram :: BinningStrategy -> [[Double]] -> Histogram
histogram strat xs = histogramNumBins (strat $ concat xs) xs

-- | Create a histogram by specifying the exact bin size. 
-- You probably don't want to use this function, and should use histogram 
-- with an appropriate binning strategy.
histogramBinSize :: Double -> [[Double]] -> Histogram
histogramBinSize size = combineBins . map (histBins size . bin size)

-- | Create a histogram by the specified approximated number of bins.
-- You probably don't want to use this function, and should use 
-- histogram with an appropriate binning strategy.
histogramNumBins :: Int -> [[Double]] -> Histogram
histogramNumBins n xs =
  histogramBinSize size xs
    where
        size = fromIntegral (firstdigit diff) * (10 ** fromIntegral (exponent10 diff))
        diff = if diff_test > 0
               then diff_test
               else 1
        diff_test = (maximum (map maximum xs) - 
                     minimum (map minimum xs)) / fromIntegral (max 1 n)

        firstdigit dbl = floor $ dbl / (10 ** fromIntegral (exponent10 dbl))
        exponent10 dbl = floor $ logBase 10 dbl

-- | It does all the binning for the histogram.
histBins :: Double -> [Double] -> [(Double, Int)]
histBins size xs = [ (head l, length l) | l <- group (sort xs) ]

-- | It "rounds" every number into the closest number below it 
-- that is divisible by size.
bin :: Double -> [Double] -> [Double]
bin size xs @ [_] = xs
bin size xs = map (\x -> size * fromIntegral (floor (x / size)) + size / 2) xs

-- | Combine bins from different histograms (optimized version).
combineBins :: [[(Double, Int)]] -> [(Double, [Int])]
combineBins [xs] = map (\(t, n) -> (t, [n])) xs
combineBins xss  = combineBins' xss

-- | Combine bins from different histograms (generic version).
combineBins' :: [[(Double, Int)]] -> [(Double, [Int])]
combineBins' xs = map (merging . sorting) $ 
                  groupping $ ordering $ concat $ numbering xs
  where numbering       = zipWith (curry indexing) [1..]
        indexing (i, x) = map (\(t, n) -> (i, t, n)) x
        ordering  = sortBy  $ \(_, t1, _) (_, t2, _) -> compare t1 t2
        groupping = groupBy $ \(_, t1, _) (_, t2, _) -> t1 == t2
        sorting   = sortBy  $ \(i1, _, _) (i2, _, _) -> compare i1 i2
        merging zs @ ((_, t, _) : _) = merging' zs t 1 []
        merging' [] t k acc
          | k <= count = merging' [] t (k + 1) (0 : acc)
          | k > count  = (t, reverse acc)
        merging' zs @ ((i, _, n) : xs) t k acc
          | i == k = merging' xs t (k + 1) (n : acc)
          | i > k  = merging' zs t (k + 1) (0 : acc)
        count = length xs

------------------------------------------------------------------------------
-- simple math functions
-- taken from package hstats, which wouldn't fully compile

-- | Numerically stable mean.
mean :: Floating a => [a] -> a
mean x = fst $ foldl' (\(m, n) x -> (m+(x-m)/(n+1),n+1)) (0,0) x

-- | Standard deviation of sample.
stddev :: (Floating a) => [a] -> a
stddev xs = sqrt $ var xs

-- | Sample variance.
var :: (Floating a) => [a] -> a
var xs = var' 0 0 0 xs / fromIntegral (length xs - 1)
    where
      var' _ _ s [] = s
      var' m n s (x:xs) = var' nm (n + 1) (s + delta * (x - nm)) xs
         where
           delta = x - m
           nm = m + delta / fromIntegral (n + 1)

-- | kurtosis is taken from wikipedia's definition.
kurtosis :: (Floating a) => [a] -> a
kurtosis xs = ((1/n) * sum [(x-x_bar)^4 | x <- xs])
            / ((1/n) * sum [(x-x_bar)^2 | x <- xs])^2 -3
    where 
        n = fromIntegral $ length xs
        x_bar = mean xs