-
Notifications
You must be signed in to change notification settings - Fork 2
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add function for generating Float/Double in [0,1] range #149
Add function for generating Float/Double in [0,1] range #149
Conversation
There're functions for Float and Double. They come in two variants: one that sample in [0,1] interval and one that samples in (0,1]. Former allows to save few floating poit operations when one needs values in [0,1] range. Which is very commonly used primitive Latter isimportant when output of generator is fed into function which not defined at 0. This is again very common. For Float sampling zero is feasible. For Double it's improbable but could in principle happen in bad generator/seed combo, or for malicious generator inicialization.
m = fromIntegral (maxBound :: Word64) :: Double | ||
{-# INLINE word64ToDoubleInUnitInterval #-} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Have you run the benchmarks on this? I have a feeling that uniformDouble01M
might need to be INLINE
to avoid a regression, but I'm not sure.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes. No change at all. I played a bit with strictness but GHC but I wasn't able to find any measurable change.
src/System/Random/Internal.hs
Outdated
|
||
-- | Generate uniformly distributed 'Double' in the range (0, 1]. That | ||
-- is result is guaranteed to be positive | ||
uniformDoublePos01M :: StatefulGen g m => g -> m Double |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
When I see "pos", I think "position". How about uniformDoublePositive01M
? Longer but unambiguous.
Additional thought: "positive 0" actually refers to +0
, as opposed to "negative zero", -0
, in IEEE 754. So people might think that something called uniformDoublePos01M
or uniformDoublePositive01M
includes positive zero. To avoid "positive" and "0" following each other in the name of the function, perhaps the scheme could be
uniformDouble01
uniformDouble01Positive
uniformFloat01
uniformFloat01Positive
-- We add small constant to shift generated value from zero. It's | ||
-- selected as 1/2 of smallest possible nonzero value | ||
d = 2.710505431213761e-20 -- 2**(-65) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is the largest value that can be added to 1.0 such that the result is still 1.0. There are many smaller values that also have this property and that would also result in never getting 0 as a result of this function. My guess is that the goal here was to avoid "excessively small values" rather than just avoiding zero, is that correct? If so, I would suggest extending the description a little to justify why this is 2^-64
and not e.g. 5e-324
, which would also avoid getting zero but would not exclude very small numbers from being returned.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It' good question.
No. It's 1/2 of smallest non zero value that this algorithm could generate. Largest ε that 1+ε=1 is 2^{-53}. Choice is somewhat arbitrary. Another possibility is to pick ε*2^{-64}. That is largest number that will leave all other values unchanged.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Right, thanks for the correction.
I guess the question is: under what circumstances do people want to avoid getting a zero here, and under those circumstances, would they also want to avoid getting extremely small values?
My default here would be to exclude only the zero. Everything else should have a good reason IMO, since excluding more numbers results in a larger effect on the outcome distribution.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
When result is fed into function that isn't defined for 0. Most common choices are: log
1/x^n
. For log anything nonzero goes. For inverse x it's better to get some distance between underflow and 1/2^64. Denormalized numbers are sort of useless here since 1/x will yield same ∞ as zero.
ε*2^{-64}. looks like good choice since it doesn't change rest of distribution and still quite far from underflows
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
1/2 of smallest possible nonzero value
I don't think that's right. The smallest positive double, or the smallest absolute value that is non-zero, is 5e-324
. So "1/2 of smallest possible nonzero value" would be 2.5e-324
.
I don't have a very strong opinion on what this number should be, but the comment should be updated.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I can explain that what @Shimuuar suggests is actually very good. Suppose you have floating point numbers with only 2 bits for the mantissa and 2 bits for the exponent then using the proposed algorithm:
word32ToFloatInUnitInterval :: Int -> Float
word32ToFloatInUnitInterval w32 = f / m
where
f = fromIntegral w32 :: Float
m = fromIntegral (maxBound :: Int) :: Float
and manually calculating for our tiny floating point values
map ((/4) . fromIntegral) [0,1,2,3]
[0.0,0.25,0.5,0.75]
The smallest value the RNG algorithm can generate is 0.25 and half of this is 0.125 so we will get
map ((+(0.25/2)) . (/4) . fromIntegral) [0,1,2,3]
[0.125,0.375,0.625,0.875]
which has the nice property that it is symmetrical in the interval. N.B. that these really are available in the floating point representation
normals :: [Double]
normals = [ 1/2^e * (1 + m0 * (1/2^1) + m1 * (1/2^2)) | e <- [2,1], m0 <- [0,1], m1 <- [0,1] ]
[0.25,0.3125,0.375,0.4375,0.5,0.625,0.75,0.875]
so no rounding is happening (that was what I was worried about).
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Of course this nice property doesn't really manifest itself with Double
and even Float
- you would have to sample a lot of Double
to notice the mean is now correct :)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ah okay that clears it up.
1/2 of smallest possible nonzero value
"smallest possible" here does not refer to the IEEE 754 encoding but to the division algorithm we use to generate doubles. Cool. Let's put that in the comment.
For concreteness' sake, here is my suggestion:
Add a small constant shift away from zero. This constant is 1/2 of the smallest positive value generated by 'uniformDouble01M'.
(And something similar for Float).
Co-authored-by: Leonhard Markert <[email protected]>
Co-authored-by: Leonhard Markert <[email protected]>
I want to check how this changes this distribution for 8 bit floating point numbers just to convince myself.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I can only recommend adding two minor things:
- A note to haddock for
uniformFloat01M
/uniformDouble01M
saying that they are used for implementinguniformM
forFloat
/Double
respectfully. @since 1.2.0
for all newly added functions in this PR
218e036
to
46b8aec
Compare
@lehins I've added since annotations and expanded documentation for Unrelated. Should we add |
Yes, we should add it for all new functions except the ones that are only exported from the |
lol. @SInCE guy probably hates by now all the Haskellers for tagging him in all the issues on github 😄 |
Instances certainly don't have such annotations. It could've been worse. Like being tagged by java programmers. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@Shimuuar I assume extra constant addition has almost no effect on performance?
I'll benchmarks for new generators and we'll compare them and uniformRM |
First round of benchmarks is very confusing. It seems we mostly measure overhead of threading state:
benchmark for Even worse. benchmark runtime jumps randomly jump when Main.hs. Few benchmarks were added and yet slow benchmak become fast and slow one fast.
It seems that we're at mercy of optimizer here. |
Also use LaTeX. It looks prettier
I've added benchmarks which thread state manually (via
What's more worrying GHC is not able to cut through state monad at -O1 and we get about 10x slowdown (haskell is usually compiled with -O1, not -O2):
|
Those are very valuable findings! Benchmarking the |
I thought |
I don't know about stack but cabal uses -O by default. It also complains about -O2. Here is what cabal thinks about math-functions:
|
Stack uses Cabal the library underneath and it does the same thing by default: |
Just to check my mental model: if I'm writing code in some project, and I depend on Related - do stack and cabal respect the Line 98 in 80a4a8a
|
It is a bad practice to add optimization flags in the library section, in fact I think Cabal will give you a warning about it if you do. It is the end executable that really matters. So if you have a project In order for inlining to be triggered |
I changed flags in the benchmark but didn't touch random library. So yes, client code should enable -O2.
It depends. Anything that is not inlined in user code will not benefit from -O2. Which may or may not be a problem P.S. I think that this PR should be merged. It didn't introduce any regressions. Only showed that we may have some performance problems unrelated to it |
[ bench "uniformFloat01M" $ nfIO $ runStateGenT_ (mkStdGen 1337) $ \g -> | ||
replicateM_ sz $ do !_ <- uniformFloat01M g |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is it possible to make this more generic, i.e. to run all benchmarks with runStateGenT_
rather than just this specific one? I would really like us to have a full benchmark matrix.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Surely it's possible. I however would prefer to do it in separate PR in order to one thing at a time
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Cool. Filed #161.
I'd still like to see #149 (comment) addressed. Otherwise LGTM! |
Fancy :) |
There're functions for Float and Double. They come in two variants: one that
sample in [0,1] interval and one that samples in (0,1]. Former allows to save
few floating poit operations when one needs values in [0,1] range. Which is very
commonly used primitive
Latter isimportant when output of generator is fed into function which not
defined at 0. This is again very common. For Float sampling zero is
feasible. For Double it's improbable but could in principle happen in bad
generator/seed combo, or for malicious generator initialization
Created as result of discussion in haskell#62