[Haskell] Random matrices

[Donald Bruce Stewart]

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

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