I'm writing some code to generate a dither (=noise) signal. I'm trying to generate an infinite series of noise with triangular distribution but my code hangs into an infinite loop. The problem is that I'm not very good with Haskell IO yet and I can't figure out how to write this piece of IO code without it looping infinitely. So, in short, how do I do this without getting into an infinite loop: tpdfs :: (Int, Int) -> IO [Int] tpdfs (low, high) = do first <- getStdRandom (randomR (low, high)) second <- getStdRandom (randomR (low, high)) let r = (first + second) `div` 2 rest <- tpdfs (low, high) return (r : rest) Caller site: do nums <- tpdfs (2, 12) let ns = take 7 nums Niko
On 7/16/07, Niko Korhonen
I'm writing some code to generate a dither (=noise) signal. I'm trying to generate an infinite series of noise with triangular distribution but my code hangs into an infinite loop. The problem is that I'm not very good with Haskell IO yet and I can't figure out how to write this piece of IO code without it looping infinitely.
So, in short, how do I do this without getting into an infinite loop:
tpdfs :: (Int, Int) -> IO [Int] tpdfs (low, high) = do first <- getStdRandom (randomR (low, high)) second <- getStdRandom (randomR (low, high)) let r = (first + second) `div` 2 rest <- tpdfs (low, high) return (r : rest)
Caller site:
do nums <- tpdfs (2, 12) let ns = take 7 nums
Niko
I did not look at it long enough to tell you why there is an infinite loop. However, think about it on a high level with me. You want a stream of these random numbers (I'm not sure what a triangular distribution is, but that's okay). To get one of these, you take two random numbers and perform a combination function (\x y -> (x + y) `div` 2 ) on them. So you can lift this from one random numbers to a stream of random numbers if you have have two streams of random numbers instead of just two random numbers. zipWith is the function that brings us from one number to a stream of numbers. tpdfs range = do g <- newStdGen -- get a random generator (g1, g2) <- return $ split g -- make two random generators out of it return $ zipWith combine (randomRs range g1) (randomRs range g2) -- get two streams of random numbers, and combine them elementwise. combine x y = (x + y) `div` 2 Uh, I know that's a very poor explanation, but hopefully it gives you an alternate way to look at the problem.
Bryan Burgers wrote:
Uh, I know that's a very poor explanation, but hopefully it gives you an alternate way to look at the problem.
Yes, this was extremely helpful, thank you very much. The moments where one realizes that a large piece of clumsy code can be replaced with a simple high-level function application seem to be an integral part of learning Haskell. This time it was zipWith. Previously (for me) it has been the folds :) I know that in Haskell there almost always is a high-level solution to a recursive problem (the legendary "here's a one-line fold that replaces your entire program"), but sometimes it can be very difficult to see, especially if IO is involved. Niko
Bryan Burgers wrote:
I did not look at it long enough to tell you why there is an infinite loop. However, think about it on a high level with me.
You want a stream of these random numbers (I'm not sure what a triangular distribution is, but that's okay). To get one of these, you take two random numbers and perform a combination function (\x y -> (x + y) `div` 2 ) on them.
Yes, precisely. Triangular distribution is a probability distribution that is equivalent to two rolls of dice. This means that the numbers at the middle of the range are much more likely to pop up than numbers at the edge of the range. It is quite close to gaussian distribution. I'm toying around with a signal processing toolkit in Haskell. The noise I'm trying to generate here is needed for a process called dithering, in which some noise is added to a quantized signal in order to improve it's accuracy. But not just any kind of noise will do for this purpose. The best noise for dithering is noise with triangular or gaussian probability distribution, instead of white noise which has equal probability distribution. But, like you said, that's not really important for the purposes of this discussion. What is is that we take a bunch of random numbers, perform some mathematical operation on them in order to introduce some statistical properties to the series and return the processed series. There are several different probability distribution functions, triangular being one of them. Triangular distribution requires two random numbers to generate one, and some functions require more than that.
So you can lift this from one random numbers to a stream of random numbers if you have have two streams of random numbers instead of just two random numbers. zipWith is the function that brings us from one number to a stream of numbers.
tpdfs range = do g <- newStdGen -- get a random generator (g1, g2) <- return $ split g -- make two random generators out of it return $ zipWith combine (randomRs range g1) (randomRs range g2) -- get two streams of random numbers, and combine them elementwise.
combine x y = (x + y) `div` 2
So, moving on to the next question, how well do you think this solution would scale if we would need n random numbers to generate one? Niko
Niko Korhonen wrote:
Bryan Burgers wrote:
tpdfs range = do g <- newStdGen -- get a random generator (g1, g2) <- return $ split g -- make two random generators out of it return $ zipWith combine (randomRs range g1) (randomRs range g2) -- get two streams of random numbers, and combine them elementwise.
combine x y = (x + y) `div` 2
So, moving on to the next question, how well do you think this solution would scale if we would need n random numbers to generate one?
As a rule of thumb, anything that works with two operands of the same type can be generalized to more than two operands, using lists and fold resp. unfold. The following code is untested and may contain errors, but you should get the idea: -- combine n generators by taking the arithmetic mean arithMean xs = sum xs `div` length xs -- must be a finite nonempty list! -- create infinitely many generators from one splitMany = map fst . iterate (split . snd) tpdfs range = do g <- newStdGen -- get a random generator gs <- return $ iterate splitMany g return $ map arithMean . map (randomRs range) (take 42 gs) Cheers Ben
Niko Korhonen wrote:
So, in short, how do I do this without getting into an infinite loop:
tpdfs :: (Int, Int) -> IO [Int] tpdfs (low, high) = do first <- getStdRandom (randomR (low, high)) second <- getStdRandom (randomR (low, high)) let r = (first + second) `div` 2 rest <- tpdfs (low, high) -- (A) return (r : rest) -- (B)
(A) will be executed before (B) because of the IO monad. But you want r to be returned before rest is computed. I would split tpdfs in two functions: a pure function converting a infinite list of random numbers to another infinite list of random numbers, and an IO-function creating the original infinite list of random numbers: tpdfs' :: [Int] -> [Int] tpdfs' (x:y:rest) = (x + y) `div` 2 : tpdfs' rest tpdfs :: (Int, Int) -> IO [Int] tpdfs range = do gen <- newStdGen return (tpdfs' (randomRs range gen)) The only IO action (newStdGen) is executed when tpdfs is called, but the infinite result list is lazily created when needed. This is possible because newStdGen uses split to create a new source of randomness exclusively for the tpdfs' function wich is not accessed anywhere else. tpdfs can be written more concisely as one of these tpdfs range = liftM (tpdfs' . randomRs range) newStdGen tpdfs range = return (tpdfs' . randomRs range) `ap` newStdGen tpdfs range = newStdGen >>= (randomRs range >>> tpdfs' >>> return) using either Control.Monad or Control.Arrow. I'm not sure your aproach is numerically correct. Let's assume range = (0, 1). The resulting number could be (0 + 0) `div` 2 = 0 (0 + 1) `div` 2 = 0 (1 + 0) `div` 2 = 0 (1 + 1) `div` 2 = 1 with equal probability. Is this what you want? Tillmann Rendel
Tillmann Rendel wrote:
(A) will be executed before (B) because of the IO monad. But you want r to be returned before rest is computed. I would split tpdfs in two functions: a pure function converting a infinite list of random numbers to another infinite list of random numbers, and an IO-function creating the original infinite list of random numbers:
tpdfs' :: [Int] -> [Int] tpdfs' (x:y:rest) = (x + y) `div` 2 : tpdfs' rest
tpdfs :: (Int, Int) -> IO [Int] tpdfs range = do gen <- newStdGen return (tpdfs' (randomRs range gen))
This seems like a good solution since it scales well with probability functions that require more than two random numbers in order to produce one. I could just write different versions of tpdfs' to process the stream and feed the functions to tpdfs. Nice.
I'm not sure your aproach is numerically correct. Let's assume range = (0, 1). The resulting number could be
(0 + 0) `div` 2 = 0 (0 + 1) `div` 2 = 0 (1 + 0) `div` 2 = 0 (1 + 1) `div` 2 = 1
with equal probability. Is this what you want?
Come to think of it, a better formula would be something like: round(x/2 + y/2) round(0/2 + 0/2) = 0 round(0/2 + 1/2) = round(0.5) = 1 round(1/2 + 0/2) = round(0.5) = 1 round(1/2 + 1/2) = = 1 But that's only because of rounding issues. Otherwise this is exactly what I want. Triangular distribution is equivalent to two rolls of dice, meaning that numbers at the middle of the range are much more likely to pop up than numbers at the edge of the range. Quite like in gaussian probability distribution. It's just that the range (0, 1) is too short for the function to work properly with integer arithmetics. It's difficult to say what the middle of the range (0, 1) should be. If we always round the result to the nearest integer, the "middle of the range" is 1. The fixed formula demostrates that the numbers at the middle of the range (round(0.5)) are most likely to appear. Niko
participants (4)
-
Benjamin Franksen -
Bryan Burgers -
Niko Korhonen -
Tillmann Rendel