[Haskell] Random matrices
Scherrer, Chad
Chad.Scherrer at pnl.gov
Sun Aug 21 02:57:09 EDT 2005
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
------------------------------------------------------------------------
More information about the Haskell
mailing list