    {-# LANGUAGE
        MultiParamTypeClasses,
        FlexibleContexts, FlexibleInstances,
        UndecidableInstances, GADTs
      #-}

    module Data.Random.Distribution.Simplex
        ( StdSimplex(..)
        , stdSimplex
        , stdSimplexT
        , fractionalStdSimplex
        ) where

    import Control.Monad
    import Data.List
    import Data.Random.RVar
    import Data.Random.Distribution
    import Data.Random.Distribution.Uniform

    -- |Uniform distribution over a standard simplex.
    newtype StdSimplex as =
        -- | @StdSimplex k@ constructs a standard simplex of dimension @k@
        -- (standard /k/-simplex).
        -- An element of the simplex represents a vector variable @as = (a_0,
        -- a_1, ..., a_k)@. The elements of @as@ are more than or equal to @0@
        -- and @sum as@ is always equal to @1@.
        StdSimplex Int
        deriving (StdSimplex as -> StdSimplex as -> Bool
    (StdSimplex as -> StdSimplex as -> Bool)
    -> (StdSimplex as -> StdSimplex as -> Bool) -> Eq (StdSimplex as)
    forall as. StdSimplex as -> StdSimplex as -> Bool
    forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
    $c== :: forall as. StdSimplex as -> StdSimplex as -> Bool
    == :: StdSimplex as -> StdSimplex as -> Bool
    $c/= :: forall as. StdSimplex as -> StdSimplex as -> Bool
    /= :: StdSimplex as -> StdSimplex as -> Bool
    Eq, Int -> StdSimplex as -> ShowS
    [StdSimplex as] -> ShowS
    StdSimplex as -> String
    (Int -> StdSimplex as -> ShowS)
    -> (StdSimplex as -> String)
    -> ([StdSimplex as] -> ShowS)
    -> Show (StdSimplex as)
    forall as. Int -> StdSimplex as -> ShowS
    forall as. [StdSimplex as] -> ShowS
    forall as. StdSimplex as -> String
    forall a.
    (Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
    $cshowsPrec :: forall as. Int -> StdSimplex as -> ShowS
    showsPrec :: Int -> StdSimplex as -> ShowS
    $cshow :: forall as. StdSimplex as -> String
    show :: StdSimplex as -> String
    $cshowList :: forall as. [StdSimplex as] -> ShowS
    showList :: [StdSimplex as] -> ShowS
    Show)

    instance (Ord a, Fractional a, Distribution StdUniform a) => Distribution StdSimplex [a] where
        rvar :: StdSimplex [a] -> RVar [a]
    rvar (StdSimplex Int
    k) = Int -> RVar [a]
    forall a.
    (Ord a, Fractional a, Distribution StdUniform a) =>
    Int -> RVar [a]
    fractionalStdSimplex Int
    k

    -- |@stdSimplex k@ returns a random variable being uniformly distributed over
    -- a standard simplex of dimension @k@.
    stdSimplex :: Distribution StdSimplex [a] => Int -> RVar [a]
    stdSimplex :: forall a. Distribution StdSimplex [a] => Int -> RVar [a]
    stdSimplex Int
    k = StdSimplex [a] -> RVar [a]
    forall (d :: * -> *) t. Distribution d t => d t -> RVar t
    rvar (Int -> StdSimplex [a]
    forall as. Int -> StdSimplex as
    StdSimplex Int
    k)

    stdSimplexT :: Distribution StdSimplex [a] => Int -> RVarT m [a]
    stdSimplexT :: forall a (m :: * -> *).
    Distribution StdSimplex [a] =>
    Int -> RVarT m [a]
    stdSimplexT Int
    k = StdSimplex [a] -> RVarT m [a]
    forall (n :: * -> *). StdSimplex [a] -> RVarT n [a]
    forall (d :: * -> *) t (n :: * -> *).
    Distribution d t =>
    d t -> RVarT n t
    rvarT (Int -> StdSimplex [a]
    forall as. Int -> StdSimplex as
    StdSimplex Int
    k)

    -- |An algorithm proposed by Rubinstein & Melamed (1998).
    -- See, /e.g./, S. Onn, I. Weissman.
    -- Generating uniform random vectors over a simplex with implications to
    -- the volume of a certain polytope and to multivariate extremes.
    -- /Ann Oper Res/ (2011) __189__:331-342.
    fractionalStdSimplex :: (Ord a, Fractional a, Distribution StdUniform a) => Int -> RVar [a]
    fractionalStdSimplex :: forall a.
    (Ord a, Fractional a, Distribution StdUniform a) =>
    Int -> RVar [a]
    fractionalStdSimplex Int
    k = do us <- Int -> RVarT Identity a -> RVarT Identity [a]
    forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
    replicateM Int
    k RVarT Identity a
    forall a. Distribution StdUniform a => RVar a
    stdUniform
                                let us' = [a] -> [a]
    forall a. Ord a => [a] -> [a]
    sort [a]
    us [a] -> [a] -> [a]
    forall a. [a] -> [a] -> [a]
    ++ [a
    1]
                                return $ zipWith (-) us' (0 : us')