I'm working with the packages System.Random.Mersenne.Pure64 and Control.Monad.Mersenne.Random by Don Stewart which are usually blazingly fast, and are supposed to help avoid common errors, like using non-strict state monads.
Nevertheless, I managed to write some code that results in a stack overflow for moderately large vectors.
import qualified Data.Vector.Unboxed as U
import Data.Int
import System.Random.Mersenne.Pure64
import Control.Monad.Mersenne.Random
main = do
    let dim = 1000000
        y = evalRandom (U.replicateM dim getInt64) (pureMT 13) :: U.Vector Int64
    putStr $ (show $ U.head y)
I'm guessing this must be due to laziness in Vector's replicateM implementation, though it's difficult to see since it is implemented using streams.
How might I write code that uses constant stack space to sample large vectors?
Nothing was immediately obviously wrong with monad-mersenne-random (didn't look at core one iota) but it's worth noting that everything works when using the state monad:
import Control.Monad.State
...
        y      = evalState (U.replicateM dim (state $ \s -> randomInt64 s)) (pureMT 13) :: U.Vector Int64
With a result of:
$ ghc --version
The Glorious Glasgow Haskell Compilation System, version 7.6.3
$ ghc so.hs -fforce-recomp
[1 of 2] Compiling Control.Monad.Mersenne.Random ( Control/Monad/Mersenne/Random.hs, Control/Monad/Mersenne/Random.o )
[2 of 2] Compiling Main             ( so.hs, so.o )
Linking so ...
$ ./so
-7188968464842378225
EDIT:
Looking at the two monad definitions (StateT and Rand) it seems the only real difference is in the strictness of the tuple.  So I tried using Control.Monad.State.Strict and ah-ha!  The stack overflow returned.  So I'm guessing that buried in Vector's replicateM code is something similar to the foldr used in the replicateM from base.  This would explain why you won't want the sequence to be strict.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With