-
-
Notifications
You must be signed in to change notification settings - Fork 433
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
Implement sampling from the unit sphere and circle #567
Conversation
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.
Nice addition.
src/distributions/unit_sphere.rs
Outdated
loop { | ||
let (x1, x2) = (self.uniform.sample(rng), self.uniform.sample(rng)); | ||
let (x1_sq, x2_sq) = (x1*x1, x2*x2); | ||
let sum = x1_sq + x2_sq; |
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.
You could just use a single sum_sq = ...
since the individual squares are not used.
src/distributions/unit_sphere.rs
Outdated
continue; | ||
} | ||
let factor = (1.0_f64 - x1_sq - x2_sq).sqrt(); | ||
return [2.*x1 * factor, 2.*x2 * factor, 1. - 2.*sum]; |
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.
pre-multiply factor by 2?
src/distributions/mod.rs
Outdated
@@ -190,6 +193,7 @@ pub use self::bernoulli::Bernoulli; | |||
pub mod uniform; | |||
mod bernoulli; | |||
#[cfg(feature="alloc")] mod weighted; | |||
#[cfg(feature="std")]mod unit_sphere; |
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.
missing a space
Thanks for the feedback, I addressed your comments. Do we want to add sampling from the unit sphere volume as well? The implementation is fairly straightforward (not very surprisingly, rejection sampling is the most efficient). If yes, I think |
To be honest I'm out of my depth regarding which sampling algorithms we should choose to implement. But I read the paper for this one and it looks like a nice approach to an apparently common problem. I don't think there's much use for sampling from the volume of a sphere (this is about the only thing I can think of), and it's not exactly difficult to implement rejection sampling, so I don't think we would want it. One can try to generalise sampling from the surface to other dimensions, but for 2D there's a choice of So the only important questions are, I think,
|
I implemented all the common ones and picked the fastest one. I think speed is the only concern, they should all produce correct uniform distributions.
FWIW, I already implemented it. I think we should include it, because the efficient implementation is non-trivial.
Sampling from the volume or the surface? I can't spontaneously think of any use cases. For large dimensions, rejection sampling becomes inefficient, so it might be useful to provide an implementation. |
Only slightly surprising. Trig is slow. Do you want to try the method I suggested? Problem is getting the sign bit without an extra RNG call requires re-implementing the conversion to FP (as in our Ziggurat function). |
It is fast, but it is unfortunately not uniform.
Maybe we should implement |
I have a uniformity test using histograms. Should I use an external crate for the histogram implementation, or should we ship our own (private) implementation in rand? |
Oh, of course not, sorry. Replacing
Not really sure what you mean here. An extra lib for testing using histograms? That sounds fairly specific to testing random number distributions so I suppose it could be a Rand sub-lib, but it could work just as well as a separate crate too. But I should probably leave it to your judgement. |
Yes, that works. It is twice as fast as the naive trigonometric approach, but still slower than the current implementation.
No, I was using an extra lib just for histograms, because I already implemented them there. In principle I could copy-and-paste the implementation into Rand's test directory. |
I don't see a problem with adding a dev-dependency on |
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.
Looks very good. The only thing I might add is that we should add a test to detect value-breaking changes. (Just recording a few outputs from the current implementations with a fixed random stream is enough.)
Is there a reason why this is always |
Perhaps I should have been clearer: we're sampling from the edge of the unit circle. (Since we return a pair of coordinates, it would be easy to mis-understand this as sampling from the disc.) In this case, it also seems there are people who want to sample from the disc itself, though I don't know if this has an application or is just a toy problem. But if we ever want to implement that, we can use the name @vks can you squash please? I don't think this needs 6 commits.
Status quo. See #100. I'd rather get this right once and break stuff as necessary than use a piecemeal solution. (Also, this won't actually be released until 0.6, so there's time to do so before it becomes a "breaking" change.) |
Oh, and CI is fixed now... but is tripping up on |
The Unity engine offers sampling from the inside of the unit circle and sphere, so I guess it is used in 3D games.
The language that I was using is the following:
I thought that would be fairly clear, but since we don't have I'll look into the remaining failure and rebasing/squashing again. |
This uses a method by Marsaglia (1972), which was found to be the fastest of the methods given by [MathWorld]. [MathWorld]: http://mathworld.wolfram.com/SpherePointPicking.html
This introduces an external dev-dependency on the `average` crate.
I clarified we sample from the edge of the circle, fixed the type inference failure and rebased/squashed. |
Great. We could rename |
This uses a method by Marsaglia (1972), which was found to be the fastest of the methods given by MathWorld. (See my benchmarks.)
There are some open questions:
UnitSphere
is slightly incorrect. We are not sampling from the the unit sphere, but rather its surface. However,UnitSphereSurface
is a bit long and sampling from the unit sphere is trivial given a sample from the surface, so I'm not sure we should add it.assert_almost_eq
could probably be used in all tests.