{-# LANGUAGE CPP #-}
{-
	Copyright (C) 2011-2017 Dr. Alistair Ward

	This program is free software: you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation, either version 3 of the License, or
	(at your option) any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with this program.  If not, see <http://www.gnu.org/licenses/>.
-}
{- |
 [@AUTHOR@]	Dr. Alistair Ward

 [@DESCRIPTION@]	Miscellaneous statistics functions.
-}

module Factory.Math.Statistics(
-- * Functions
        getMean,
        getRootMeanSquare,
        getWeightedMean,
--	getDispersionFromMean,
        getVariance,
        getStandardDeviation,
        getAverageAbsoluteDeviation,
        getCoefficientOfVariance,
        nCr,
        nPr
) where

import                  Control.Arrow((***))
import qualified        Control.Exception
import                  Control.Parallel(par, pseq)
import qualified        Data.Foldable
import qualified        Data.List
import qualified        Factory.Math.Factorial                  as Math.Factorial
import qualified        Factory.Math.Implementations.Factorial  as Math.Implementations.Factorial
import qualified        Factory.Math.Power                      as Math.Power

{- |
	* Determines the /mean/ of the specified numbers; <https://en.wikipedia.org/wiki/Mean>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getMean :: (
        Data.Foldable.Foldable  foldable,
        Fractional              result,
        Real                    value
 )
        => foldable value
        -> result
{-# SPECIALISE getMean :: Data.Foldable.Foldable foldable => foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getMean :: Data.Foldable.Foldable foldable => foldable Rational -> Rational #-}
#endif
getMean foldable        = Control.Exception.assert (denominator /= 0) $ realToFrac numerator / fromIntegral denominator where
        denominator :: Int
        (numerator, denominator)        = Data.Foldable.foldl' (
                \acc x  -> let
                        acc'@(n, d)     = (+ x) *** succ $ acc
                in n `seq` d `seq` acc'
         ) (0, 0) foldable

-- | Determines the /root mean square/ of the specified numbers; <https://en.wikipedia.org/wiki/Root_mean_square>.
getRootMeanSquare :: (
        Data.Foldable.Foldable  foldable,
        Floating                result,
        Real                    value
 )
        => foldable value
        -> result
{-# SPECIALISE getRootMeanSquare :: Data.Foldable.Foldable foldable => foldable Double -> Double #-}
getRootMeanSquare foldable      = Control.Exception.assert (denominator /= 0) . sqrt $ realToFrac numerator / fromIntegral denominator  where
        denominator :: Int
        (numerator, denominator)        = Data.Foldable.foldl' (
                \acc x -> let
                        acc'@(n, d)     = (+ Math.Power.square x) *** succ $ acc
                 in n `seq` d `seq` acc'
         ) (0, 0) foldable

{- |
	* Determines the /weighted mean/ of the specified numbers; <https://en.wikipedia.org/wiki/Weighted_arithmetic_mean>.

	* The specified value is only evaluated if the corresponding weight is non-zero.

	* CAVEAT: because the operand is more general than a list, no optimisation is performed when supplied a singleton.
-}
getWeightedMean :: (
        Data.Foldable.Foldable  foldable,
        Fractional              result,
        Real                    value,
        Real                    weight
 )
        => foldable (value, weight)     -- ^ Each pair consists of a value & the corresponding weight.
        -> result
{-# SPECIALISE getWeightedMean :: Data.Foldable.Foldable foldable => foldable (Double, Double) -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getWeightedMean :: Data.Foldable.Foldable foldable => foldable (Rational, Rational) -> Rational #-}
#endif
getWeightedMean foldable        = Control.Exception.assert (denominator /= 0) $ numerator / realToFrac denominator      where
        (numerator, denominator)        = Data.Foldable.foldl' (
                \acc (value, weight)    -> if weight == 0
                        then acc        -- Avoid unnecessary evaluation.
                        else let
                                acc'@(n, d)     = (+ realToFrac weight * realToFrac value) *** (+ weight) $ acc
                        in n `seq` d `seq` acc'
         ) (0, 0) foldable

{- |
	* Measures the /dispersion/ of a /population/ of results from the /mean/ value; <https://en.wikipedia.org/wiki/Statistical_dispersion>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getDispersionFromMean :: (
        Data.Foldable.Foldable  foldable,
        Fractional              result,
        Functor                 foldable,
        Real                    value
 ) => (Rational -> Rational) -> foldable value -> result
{-# SPECIALISE getDispersionFromMean :: (Data.Foldable.Foldable foldable, Functor foldable) => (Rational -> Rational) -> foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getDispersionFromMean :: (Data.Foldable.Foldable foldable, Functor foldable) => (Rational -> Rational) -> foldable Rational -> Rational #-}
#endif
getDispersionFromMean weight foldable   = getMean $ fmap (weight . subtract mean . toRational) foldable where
        mean    = getMean foldable

{- |
	* Determines the exact /variance/ of the specified numbers; <https://en.wikipedia.org/wiki/Variance>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getVariance :: (
        Data.Foldable.Foldable  foldable,
        Fractional              variance,
        Functor                 foldable,
        Real                    value
 ) => foldable value -> variance
{-# SPECIALISE getVariance :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getVariance :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Rational -> Rational #-}
#endif
getVariance     = getDispersionFromMean Math.Power.square

-- | Determines the /standard-deviation/ of the specified numbers; <https://en.wikipedia.org/wiki/Standard_deviation>.
getStandardDeviation :: (
        Data.Foldable.Foldable  foldable,
        Floating                result,
        Functor                 foldable,
        Real                    value
 ) => foldable value -> result
{-# SPECIALISE getStandardDeviation :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
getStandardDeviation    = sqrt . getVariance

{- |
	* Determines the /average absolute deviation/ of the specified numbers; <https://en.wikipedia.org/wiki/Absolute_deviation#Average_absolute_deviation>.

	* Should the caller define the result-type as 'Rational', then it will be free from rounding-errors.
-}
getAverageAbsoluteDeviation :: (
        Data.Foldable.Foldable  foldable,
        Fractional              result,
        Functor                 foldable,
        Real                    value
 ) => foldable value -> result
{-# SPECIALISE getAverageAbsoluteDeviation :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
#if MIN_TOOL_VERSION_ghc(7,10,0)
{-# SPECIALISE getAverageAbsoluteDeviation :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Rational -> Rational #-}
#endif
getAverageAbsoluteDeviation     = getDispersionFromMean abs

-- | Determines the /coefficient-of-variance/ of the specified numbers; <https://en.wikipedia.org/wiki/Coefficient_of_variation>.
getCoefficientOfVariance :: (
        Data.Foldable.Foldable  foldable,
        Eq                      result,
        Floating                result,
        Functor                 foldable,
        Real                    value
 ) => foldable value -> result
{-# SPECIALISE getCoefficientOfVariance :: (Data.Foldable.Foldable foldable, Functor foldable) => foldable Double -> Double #-}
getCoefficientOfVariance l      = Control.Exception.assert (mean /= 0) $ getStandardDeviation l / abs mean      where
        mean    = getMean l

-- | The number of unordered /combinations/ of /r/ objects taken from /n/; <https://en.wikipedia.org/wiki/Combination>.
nCr :: (Math.Factorial.Algorithmic factorialAlgorithm, Integral i, Show i)
        => factorialAlgorithm
        -> i    -- ^ The total number of items from which to select.
        -> i    -- ^ The number of items in a sample.
        -> i    -- ^ The number of combinations.
nCr _ 0 _       = 1
nCr _ _ 0       = 1
nCr factorialAlgorithm n r
        | n < r         = 0
        | otherwise     = Control.Exception.assert (n >= 0 && r >= 0) $ numerator `par` (denominator `pseq` numerator `div` denominator)
        where
                [smaller, bigger]       = Data.List.sort [r, n - r]
                numerator               = Math.Implementations.Factorial.risingFactorial (succ bigger) (n - bigger)
                denominator             = Math.Factorial.factorial factorialAlgorithm smaller

-- | The number of /permutations/ of /r/ objects taken from /n/; <https://en.wikipedia.org/wiki/Permutations>.
nPr :: (Integral i, Show i)
        => i    -- ^ The total number of items from which to select.
        -> i    -- ^ The number of items in a sample.
        -> i    -- ^ The number of permutations.
nPr 0 _ = 1
nPr _ 0 = 1
nPr n r
        | n < r         = 0
        | otherwise     = Control.Exception.assert (n >= 0 && r >= 0) $ Math.Implementations.Factorial.fallingFactorial n r