{- | Some extra facilities to work with 'Rand' monad and 'PureMT'
     random number generator.
-}

module Moo.GeneticAlgorithm.Random
    (
    -- * Random numbers from given range

      getRandomR
    , getRandom
    -- * Probability distributions

    , getNormal2
    , getNormal
    -- * Random samples and shuffles

    , randomSample
    , randomSampleIndices
    , shuffle
    -- * Building blocks

    , withProbability
    -- * Re-exports from random number generator packages

    , getBool, getInt, getWord, getInt64, getWord64, getDouble
    , runRand, evalRand, newPureMT, liftRand
    , Rand, Random, PureMT
    ) where

import Control.Monad (liftM)
import qualified Control.Monad.Random.Strict as MonadRandom
import Control.Monad.Random.Strict (liftRand, runRand, evalRand)
import Data.Complex (Complex (..))
import Data.Int (Int64)
import Data.Word (Word64)
import System.Random (RandomGen, Random(..))
import System.Random.Mersenne.Pure64
import qualified System.Random.Shuffle as S
import qualified Data.Set as Set

type Rand = MonadRandom.Rand PureMT

-- | Yield a new randomly selected value of type @a@ in the range @(lo, hi)@.

-- See 'System.Random.randomR' for details.

getRandomR :: Random a => (a, a) -> Rand a
getRandomR :: (a, a) -> Rand a
getRandomR (a, a)
range = (PureMT -> (a, PureMT)) -> Rand a
forall g a. (g -> (a, g)) -> Rand g a
liftRand ((PureMT -> (a, PureMT)) -> Rand a)
-> (PureMT -> (a, PureMT)) -> Rand a
forall a b. (a -> b) -> a -> b
$ \PureMT
s -> (a, a) -> PureMT -> (a, PureMT)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (a, a)
range PureMT
s

-- | Yield a new randomly selected value of type @a@.

-- See 'System.Random.random' for details.

getRandom :: Random a => Rand a
getRandom :: Rand a
getRandom = (PureMT -> (a, PureMT)) -> Rand a
forall g a. (g -> (a, g)) -> Rand g a
liftRand PureMT -> (a, PureMT)
forall a g. (Random a, RandomGen g) => g -> (a, g)
random

getBool :: Rand Bool
getBool :: Rand Bool
getBool = Rand Bool
forall a. Random a => Rand a
getRandom
getDouble :: Rand Double
getDouble :: Rand Double
getDouble = Rand Double
forall a. Random a => Rand a
getRandom
getWord :: Rand Word
getWord :: Rand Word
getWord = Rand Word
forall a. Random a => Rand a
getRandom
getInt :: Rand Int
getInt :: Rand Int
getInt = Rand Int
forall a. Random a => Rand a
getRandom
getInt64 :: Rand Int64
getInt64 :: Rand Int64
getInt64 = Rand Int64
forall a. Random a => Rand a
getRandom
getWord64 :: Rand Word64
getWord64 :: Rand Word64
getWord64 = Rand Word64
forall a. Random a => Rand a
getRandom

-- | Yield two randomly selected values which follow standard

-- normal distribution.

getNormal2 :: Rand (Double, Double)
getNormal2 :: Rand (Double, Double)
getNormal2 = do
  -- Box-Muller method

  Double
u <- Rand Double
getDouble
  Double
v <- Rand Double
getDouble
  let (Double
c :+ Double
s) = Complex Double -> Complex Double
forall a. Floating a => a -> a
exp (Double
0 Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
v))
  let r :: Double
r = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (-Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
u
  (Double, Double) -> Rand (Double, Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
c, Double
rDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
s)

-- | Yield one randomly selected value from standard normal distribution.

getNormal :: Rand Double
getNormal :: Rand Double
getNormal = (Double, Double) -> Double
forall a b. (a, b) -> a
fst ((Double, Double) -> Double)
-> Rand (Double, Double) -> Rand Double
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
`liftM` Rand (Double, Double)
getNormal2

-- | Take at most n random elements from the list. Preserve order.

randomSample :: Int -> [a] -> Rand [a]
randomSample :: Int -> [a] -> Rand [a]
randomSample Int
n [a]
xs =
  (PureMT -> ([a], PureMT)) -> Rand [a]
forall g a. (g -> (a, g)) -> Rand g a
liftRand ((PureMT -> ([a], PureMT)) -> Rand [a])
-> (PureMT -> ([a], PureMT)) -> Rand [a]
forall a b. (a -> b) -> a -> b
$ \PureMT
g -> PureMT -> Int -> Int -> [a] -> [a] -> ([a], PureMT)
forall t a.
RandomGen t =>
t -> Int -> Int -> [a] -> [a] -> ([a], t)
select PureMT
g Int
n ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) [a]
xs []
  where
    select :: t -> Int -> Int -> [a] -> [a] -> ([a], t)
select t
rng Int
_ Int
_ [] [a]
acc = ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
acc, t
rng)
    select t
rng Int
n Int
m [a]
xs [a]
acc
        | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0     = ([a] -> [a]
forall a. [a] -> [a]
reverse [a]
acc, t
rng)
        | Bool
otherwise  =
            let (Int
k, t
rng') = (Int, Int) -> t -> (Int, t)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Int
0, Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n) t
rng
                (a
x:[a]
rest) = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
k [a]
xs
            in  t -> Int -> Int -> [a] -> [a] -> ([a], t)
select t
rng' (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) [a]
rest (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
acc)

-- | Select @sampleSize@ numbers in the range from @0@ to @(populationSize-1)@.

-- The function works best when @sampleSize@ is much smaller than @populationSize@.

randomSampleIndices :: Int -> Int -> Rand [Int]
randomSampleIndices :: Int -> Int -> Rand [Int]
randomSampleIndices Int
sampleSize Int
populationSize =
    (PureMT -> ([Int], PureMT)) -> Rand [Int]
forall g a. (g -> (a, g)) -> Rand g a
liftRand ((PureMT -> ([Int], PureMT)) -> Rand [Int])
-> (PureMT -> ([Int], PureMT)) -> Rand [Int]
forall a b. (a -> b) -> a -> b
$ \PureMT
g ->
        let (Set Int
sampleSet, PureMT
g') = PureMT -> Int -> Set Int -> (Set Int, PureMT)
forall a t.
(Ord a, Num a, RandomGen t) =>
t -> a -> Set Int -> (Set Int, t)
buildSampleSet PureMT
g Int
sampleSize Set Int
forall a. Set a
Set.empty
        in  (Set Int -> [Int]
forall a. Set a -> [a]
Set.toList Set Int
sampleSet, PureMT
g')
  where
    buildSampleSet :: t -> a -> Set Int -> (Set Int, t)
buildSampleSet t
g a
n Set Int
s
        | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
0 = (Set Int
s, t
g)
        | Bool
otherwise =
            let (Int
i, t
g') = (Int, Int) -> t -> (Int, t)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Int
0, Int
populationSizeInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) t
g
            in  if (Int
i Int -> Set Int -> Bool
forall a. Ord a => a -> Set a -> Bool
`Set.member` Set Int
s)
                then t -> a -> Set Int -> (Set Int, t)
buildSampleSet t
g' a
n Set Int
s
                else t -> a -> Set Int -> (Set Int, t)
buildSampleSet t
g' (a
na -> a -> a
forall a. Num a => a -> a -> a
-a
1) (Int -> Set Int -> Set Int
forall a. Ord a => a -> Set a -> Set a
Set.insert Int
i Set Int
s)

-- | Randomly reorder the list.

shuffle :: [a] -> Rand [a]
shuffle :: [a] -> Rand [a]
shuffle [a]
xs = (PureMT -> ([a], PureMT)) -> Rand [a]
forall g a. (g -> (a, g)) -> Rand g a
liftRand ((PureMT -> ([a], PureMT)) -> Rand [a])
-> (PureMT -> ([a], PureMT)) -> Rand [a]
forall a b. (a -> b) -> a -> b
$ \PureMT
g -> [a] -> Int -> PureMT -> ([a], PureMT)
forall gen a. RandomGen gen => [a] -> Int -> gen -> ([a], gen)
randomShuffle [a]
xs ([a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
xs) PureMT
g

-- | Given a sequence (e1,...en) to shuffle, its length, and a random

-- generator, compute the corresponding permutation of the input

-- sequence, return the permutation and the new state of the

-- random generator.

randomShuffle :: RandomGen gen => [a] -> Int -> gen -> ([a], gen)
randomShuffle :: [a] -> Int -> gen -> ([a], gen)
randomShuffle [a]
elements Int
len gen
g =
    let ([Int]
rs, gen
g') = Int -> gen -> ([Int], gen)
forall gen. RandomGen gen => Int -> gen -> ([Int], gen)
rseq Int
len gen
g
    in  ([a] -> [Int] -> [a]
forall a. [a] -> [Int] -> [a]
S.shuffle [a]
elements [Int]
rs, gen
g')
  where
  -- | The sequence (r1,...r[n-1]) of numbers such that r[i] is an

  -- independent sample from a uniform random distribution

  -- [0..n-i]

  rseq :: RandomGen gen => Int -> gen -> ([Int], gen)
  rseq :: Int -> gen -> ([Int], gen)
rseq Int
n gen
g = ([gen] -> gen) -> ([Int], [gen]) -> ([Int], gen)
forall b c a. (b -> c) -> (a, b) -> (a, c)
second [gen] -> gen
lastGen (([Int], [gen]) -> ([Int], gen))
-> ([(Int, gen)] -> ([Int], [gen])) -> [(Int, gen)] -> ([Int], gen)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Int, gen)] -> ([Int], [gen])
forall a b. [(a, b)] -> ([a], [b])
unzip ([(Int, gen)] -> ([Int], gen)) -> [(Int, gen)] -> ([Int], gen)
forall a b. (a -> b) -> a -> b
$ Int -> gen -> [(Int, gen)]
forall gen. RandomGen gen => Int -> gen -> [(Int, gen)]
rseq' (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) gen
g
      where
        rseq' :: RandomGen gen => Int -> gen -> [(Int, gen)]
        rseq' :: Int -> gen -> [(Int, gen)]
rseq' Int
i gen
gen
          | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0    = []
          | Bool
otherwise = let (Int
j, gen
gen') = (Int, Int) -> gen -> (Int, gen)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Int
0, Int
i) gen
gen
                        in  (Int
j, gen
gen') (Int, gen) -> [(Int, gen)] -> [(Int, gen)]
forall a. a -> [a] -> [a]
: Int -> gen -> [(Int, gen)]
forall gen. RandomGen gen => Int -> gen -> [(Int, gen)]
rseq' (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) gen
gen'
        -- apply a function on the second element of a pair

        second :: (b -> c) -> (a, b) -> (a, c)
        second :: (b -> c) -> (a, b) -> (a, c)
second b -> c
f (a
x,b
y) = (a
x, b -> c
f b
y)
        -- the last returned random number generator

        lastGen :: [gen] -> gen
lastGen [] = gen
g   -- didn't use the generator yet

        lastGen (gen
lst:[]) = gen
lst
        lastGen [gen]
gens = [gen] -> gen
lastGen (Int -> [gen] -> [gen]
forall a. Int -> [a] -> [a]
drop Int
1 [gen]
gens)

-- |Modify value with probability @p@. Return the unchanged value with probability @1-p@.

withProbability :: Double -> (a -> Rand a) -> (a -> Rand a)
withProbability :: Double -> (a -> Rand a) -> a -> Rand a
withProbability Double
p a -> Rand a
modify a
x = do
  Double
t <- Rand Double
getDouble
  if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
p
     then a -> Rand a
modify a
x
     else a -> Rand a
forall (m :: * -> *) a. Monad m => a -> m a
return a
x