[Haskell] Random matrices

[Scherrer, Chad]

Thanks, Don. I need some time to take a closer look at this. I'm trying to see, what is it about the code you posted that makes it so much more efficient? Is it the loop in the IO monad is compiled into an honest-to-goodness loop? Is this specific to IO, or does it work out this way for monads in general? Or maybe some gains were made by switching from using StdGen with a specified seed to using the OS-supplied random number generator?

-Chad


-----Original Message-----
From: Donald Bruce Stewart [mailto:dons at cse.unsw.edu.au]
Sent: Fri 8/19/2005 7:42 PM
To: Scherrer, Chad
Cc: haskell at haskell.org
Subject: Re: [Haskell] Random matrices
 
Chad.Scherrer:
> 
>    I'm doing some statistical calculations, and I decided to
>    try this out in Haskell to see how it goes. I'm really
>    enjoying using the language, so I hope I can straighten this
>    out so I can keep using it at work.
> 
>    I keep getting stack overflows, so I think my code must be
>    too lazy somewhere (that's what that means, right?) Here is
>    my code to build random vectors and matrices.

I was suspicious about all those zips and repeats, so did a bit of
profiling. This only revealed that the creation of random floats was
chewing up about 70% of time and space. Still suspicious of that
replicate code, I decided to rewrite it as a loop (I moved the code
into the IO monad just because it is easier to hack in IO). This fixed the
issue.

-- Don

Before:

    paprika$ ghc -package mtl -O -prof -auto-all N.hs
    paprika$ ./a.out +RTS -p
    Stack space overflow: current size 8388608 bytes.
    Use `+RTS -Ksize' to increase it.

            total time  =        0.16 secs   (8 ticks @ 20 ms)
            total alloc =  42,614,984 bytes  (excludes profiling overheads)

    COST CENTRE                    MODULE               %time %alloc
    randomEntry                    Main                  50.0   78.0
    randomArray                    Main                  50.0   19.6

------------------------------------------------------------------------

After:
    paprika$ ghc -O -prof -auto-all M.hs 
    paprika$ ./a.out +RTS -p
    array (1,10) [(1,0.73394084),(2,0.39095977),(3,0.18828022),(4,0.19094983),(5,0.83119744),(6,0.3594179),(7,0.43519533),(8,0.39867708),(9,0.15676379),(10,0.4503187)]

            total time  =        0.00 secs   (0 ticks @ 20 ms)
            total alloc =     141,368 bytes  (excludes profiling overheads)

    COST CENTRE                    MODULE               %time %alloc
    CAF                            System.Random          0.0    5.6
    CAF                            GHC.Handle             0.0    1.0
    CAF                            GHC.Float              0.0   35.2
    randomEntry                    Main                   0.0   25.7

------------------------------------------------------------------------

    import Data.Array.Unboxed
    import Data.Array.Diff
    import System.Random
    import Control.Monad

    type Vector = UArray Int Float
    type Matrix = DiffArray Int Vector

    -- Generate a random number from the unit interval
    randomEntry :: IO Float
    randomEntry = getStdRandom (randomR (0,1))

    -- Build an array of n random things
    randomArray :: (IArray arr a) => Int -> IO a -> IO (arr Int a)
    randomArray n x = do
            let loop i | i == 0    = return []
                       | otherwise = do f  <- x
                                        fs <- loop (i-1)
                                        return (f : fs)
            fs <- loop n
            return $ listArray (1,n) fs

    -- A random vector is an array of random entries
    randomVector :: Int -> IO Vector
    randomVector n = randomArray n randomEntry

    -- a random matrix is an array of random vectors.
    randomMatrix :: Int -> Int -> IO Matrix
    randomMatrix i j = randomArray i (randomVector j)

    tester :: Int -> IO [Matrix]
    tester n = liftM repeat $ randomMatrix n n

    main :: IO ()
    main = tester 10 >>= \as -> print $ (as !! 10000) ! 5

------------------------------------------------------------------------