{-
	Copyright (C) 2011 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 functions for probability-distributions.
-}

module Factory.Math.Probability(
-- * Types
-- ** Data-types
	ContinuousDistribution(..),
	DiscreteDistribution(..),
-- * Functions
	boxMullerTransform,
--	minPositiveFloat,
--	reProfile,
	generateStandardizedNormalDistribution,
	generateContinuousPopulation,
	generatePoissonDistribution,
	generateDiscretePopulation
) where

import qualified	Control.Arrow
import			Control.Arrow((***), (&&&))
import qualified	Factory.Data.Interval	as Data.Interval
import qualified	System.Random
import qualified	ToolShed.Data.List
import qualified	ToolShed.Data.Pair
import qualified	ToolShed.SelfValidate

-- | Describes a /continuous probability-distribution/; <http://en.wikipedia.org/wiki/List_of_probability_distributions#Continuous_distributions>.
data ContinuousDistribution f
	= UniformDistribution (Data.Interval.Interval f)	-- ^ Defines a /Uniform/-distribution within a closed /interval/; <http://en.wikipedia.org/wiki/Uniform_distribution>.
	| NormalDistribution f f				-- ^ Defines a /Normal/-distribution with a particular /mean/ and /variance/; <http://en.wikipedia.org/wiki/Normal_distribution>.
	deriving (Eq, Read, Show)

instance (Num a, Ord a, Show a) => ToolShed.SelfValidate.SelfValidator (ContinuousDistribution a)	where
	getErrors distribution	= ToolShed.SelfValidate.extractErrors $ case distribution of
		UniformDistribution interval	-> [(Data.Interval.isReversed interval, "Reversed interval='" ++ show interval ++ "'.")]
		NormalDistribution _ v		-> [(v < 0, "Negative variance=" ++ show v ++ ".")]

-- | Describes a /discrete probability-distribution/; <http://en.wikipedia.org/wiki/List_of_probability_distributions#Discrete_distributions>.
data DiscreteDistribution f	= PoissonDistribution f	deriving (Eq, Read, Show)

instance (Num f, Ord f, Show f) => ToolShed.SelfValidate.SelfValidator (DiscreteDistribution f)	where
	getErrors (PoissonDistribution lambda)	= ToolShed.SelfValidate.extractErrors [(lambda < 0, "Negative lambda=" ++ show lambda ++ ".")]

{- |
	* Converts a pair of independent /uniformly distributed/ random numbers, within the /semi-closed/ /unit interval/ /(0 .. 1]/,
	to a pair of independent /normally distributed/ random numbers, of standardized /mean/=0, and /variance/=1.

	* <http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform>.
-}
boxMullerTransform :: (Floating f, Ord f, Show f)
	=> (f, f)	-- ^ Independent, /uniformly distributed/ random numbers, which must be within the /semi-closed unit interval/, /(0, 1]/.
	-> (f, f)	-- ^ Independent, /normally distributed/ random numbers, with standardized /mean/=0 and /variance/=1.
boxMullerTransform cartesian
	| not . uncurry (&&) $ ToolShed.Data.Pair.mirror inSemiClosedUnitInterval cartesian	= error $ "Factory.Math.Probability.boxMullerTransform:\tspecified Cartesian coordinates, must be within semi-closed unit-interval (0, 1]; " ++ show cartesian
	| otherwise								= polarToCartesianTransform $ (sqrt . negate . (* 2) . log *** (* 2) . (* pi)) cartesian
	where
		inSemiClosedUnitInterval :: (Num n, Ord n) => n -> Bool
		inSemiClosedUnitInterval	= uncurry (&&) . ((> 0) &&& (<= 1))

		polarToCartesianTransform :: Floating f => (f, f) -> (f, f)
		polarToCartesianTransform	= uncurry (*) . Control.Arrow.second cos &&& uncurry (*) . Control.Arrow.second sin

{- |
	* Determines the minimum positive floating-point number, which can be represented by using the parameter's type.

	* Only the type of the parameter is relevant, not its value.
-}
minPositiveFloat :: RealFloat a => a -> a
minPositiveFloat	= encodeFloat 1 . uncurry (-) . (fst . floatRange &&& floatDigits)

{- |
	* Uses the supplied random-number generator,
	to generate a conceptually infinite list, of /normally distributed/ random numbers, with standardized /mean/=0, and /variance/=1.

	* <http://en.wikipedia.org/wiki/Normal_distribution>, <http://mathworld.wolfram.com/NormalDistribution.html>.
-}
generateStandardizedNormalDistribution :: (
	RealFloat		f,
	Show			f,
	System.Random.Random	f,
	System.Random.RandomGen	randomGen
 ) => randomGen -> [f]
generateStandardizedNormalDistribution	= ToolShed.Data.List.linearise . uncurry (zipWith $ curry boxMullerTransform) . ToolShed.Data.Pair.mirror (
	System.Random.randomRs (minPositiveFloat undefined, 1)
 ) . System.Random.split

-- | Stretches and shifts a /standardized normal distribution/ to achieve the required /mean/ and /standard-deviation/.
reProfile :: Num n => n -> n -> [n] -> [n]
reProfile mean standardDeviation = map ((+ mean) . (* standardDeviation))

{- |
	* Generates a random sample-population, with the specified continuous probability-distribution.

	* When a /Normal distribution/ is requested,
	the generated population will only tend towards the requested /mean/ and /variance/ of, as the sample-size tends towards infinity.
	Whilst one could arrange for these criteria to be precisely met for any sample-size, the sample would lose a degree of randomness as a result.
-}
generateContinuousPopulation :: (
	RealFloat		f,
	Show			f,
	System.Random.Random	f,
	System.Random.RandomGen	randomGen
 )
	=> Int	-- ^ number of items.
	-> ContinuousDistribution f
	-> randomGen	-- ^ A generator of /uniformly distributed/ random numbers.
	-> [f]
generateContinuousPopulation 0 _ _				= []
generateContinuousPopulation populationSize probabilityDistribution randomGen
	| populationSize < 0						= error $ "Factory.Math.Probability.generateDiscretePopulation:\tinvalid population-size=" ++ show populationSize
	| not $ ToolShed.SelfValidate.isValid probabilityDistribution	= error $ "Factory.Math.Probability.generateContinuousPopulation:\t" ++ ToolShed.SelfValidate.getFirstError probabilityDistribution
	| otherwise						= take populationSize $ (
		case probabilityDistribution of
			UniformDistribution interval				-> System.Random.randomRs interval
			NormalDistribution requiredMean requiredVariance	-> reProfile requiredMean (sqrt requiredVariance) . generateStandardizedNormalDistribution
	) randomGen

{- |
	* Uses the supplied random-number generator,
	to generate a conceptually infinite list, of random integers conforming to the /Poisson distribution/ (/mean/=lambda, /variance/=lambda).

	* <http://en.wikipedia.org/wiki/Poisson_distribution>.

	* CAVEAT:
		uses an algorithm by Knuth, which having a /linear time-complexity/ in /lambda/, can be intolerably slow;
		also, the term @exp $ negate lambda@, underflows for large /lambda/;
		so for large /lambda/, this implementation returns the appropriate 'NormalDistribution', which is similar for large /lambda/.
-}
generatePoissonDistribution :: (
	Integral		events,
	RealFloat		lambda,
	Show			lambda,
	System.Random.Random	lambda,
	System.Random.RandomGen	randomGen
 )
	=> lambda	-- ^ Defines the required approximate value of both /mean/ and /variance/.
	-> randomGen
	-> [events]
generatePoissonDistribution lambda
	| lambda < 0	= error $ "Factory.Math.Probability.generatePoissonDistribution:\tinvalid lambda=" ++ show lambda
	| lambda > (
		negate . log $ minPositiveFloat lambda	--Guard against underflow, in the user-defined type for lambda.
	)		= filter (>= 0) . map round . reProfile lambda (sqrt lambda) . generateStandardizedNormalDistribution
	| otherwise	= generator
	where
		generator	= uncurry (:) . (
			fst . head . dropWhile (
				(> exp (negate lambda)) . snd	--CAVEAT: underflows if lambda > (103 :: Float, 745 :: Double).
			) . scanl (
				\accumulator random	-> succ *** (* random) $ accumulator
			) (negate 1, 1) . System.Random.randomRs (0, 1) *** generator {-recurse-}
		 ) . System.Random.split

-- | Generates a random sample-population, with the specified discrete probability-distribution.
generateDiscretePopulation :: (
	Ord			f,
	RealFloat		f,
	Show			f,
	System.Random.Random	f,
	System.Random.RandomGen	randomGen,
	Integral		events
 )
	=> Int	-- ^ number of items.
	-> DiscreteDistribution f
	-> randomGen	-- ^ A generator of /uniformly distributed/ random numbers.
	-> [events]
generateDiscretePopulation 0 _ _				= []
generateDiscretePopulation populationSize probabilityDistribution randomGen
	| populationSize < 0						= error $ "Factory.Math.Probability.generateDiscretePopulation:\tinvalid populationSize=" ++ show populationSize
	| not $ ToolShed.SelfValidate.isValid probabilityDistribution	= error $ "Factory.Math.Probability.generateDiscretePopulation:\t" ++ ToolShed.SelfValidate.getFirstError probabilityDistribution
	| otherwise							= take populationSize $ (
		case probabilityDistribution of
			PoissonDistribution lambda	-> generatePoissonDistribution lambda
	) randomGen