Merge pull request rust-random#798 from vks/unit-ball

Implement unit ball sampling and rename UnitSphereSurface -> UnitSphere

Author: dhardy

Cargo.toml

@@ -20,3 +20,7 @@ appveyor = { repository = "rust-random/rand" }
 
 [dependencies]
 rand = { path = "..", version = ">=0.5, <=0.7" }
+
+[dev-dependencies]
+# Histogram implementation for testing uniformity
+average = "0.9.2"

src/lib.rs

@@ -55,14 +55,18 @@
 //!   - [`Triangular`] distribution
 //! - Multivariate probability distributions
 //!   - [`Dirichlet`] distribution
-//!   - [`UnitSphereSurface`] distribution
+//!   - [`UnitSphere`] distribution
+//!   - [`UnitBall`] distribution
 //!   - [`UnitCircle`] distribution
+//!   - [`UnitDisc`] distribution
 
 pub use rand::distributions::{Distribution, DistIter, Standard,
     Alphanumeric, Uniform, OpenClosed01, Open01, Bernoulli, uniform, weighted};
 
-pub use self::unit_sphere::UnitSphereSurface;
+pub use self::unit_sphere::UnitSphere;
+pub use self::unit_ball::UnitBall;
 pub use self::unit_circle::UnitCircle;
+pub use self::unit_disc::UnitDisc;
 pub use self::gamma::{Gamma, Error as GammaError, ChiSquared, ChiSquaredError,
     FisherF, FisherFError, StudentT, Beta, BetaError};
 pub use self::normal::{Normal, Error as NormalError, LogNormal, StandardNormal};
@@ -78,7 +82,9 @@ pub use self::weibull::{Weibull, Error as WeibullError};
 pub use self::utils::Float;
 
 mod unit_sphere;
+mod unit_ball;
 mod unit_circle;
+mod unit_disc;
 mod gamma;
 mod normal;
 mod exponential;

src/unit_ball.rs

@@ -0,0 +1,69 @@
+// Copyright 2019 Developers of the Rand project.
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
+// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
+// option. This file may not be copied, modified, or distributed
+// except according to those terms.
+
+use rand::Rng;
+use crate::{Distribution, Uniform, uniform::SampleUniform};
+use crate::utils::Float;
+
+/// Samples uniformly from the unit ball (surface and interior) in three
+/// dimensions.
+///
+/// Implemented via rejection sampling.
+///
+///
+/// # Example
+///
+/// ```
+/// use rand_distr::{UnitBall, Distribution};
+///
+/// let v: [f64; 3] = UnitBall.sample(&mut rand::thread_rng());
+/// println!("{:?} is from the unit ball.", v)
+/// ```
+#[derive(Clone, Copy, Debug)]
+pub struct UnitBall;
+
+impl<N: Float + SampleUniform> Distribution<[N; 3]> for UnitBall {
+    #[inline]
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [N; 3] {
+        let uniform = Uniform::new(N::from(-1.), N::from(1.));
+        let mut x1;
+        let mut x2;
+        let mut x3;
+        loop {
+            x1 = uniform.sample(rng);
+            x2 = uniform.sample(rng);
+            x3 = uniform.sample(rng);
+            if x1*x1 + x2*x2 + x3*x3 <= N::from(1.) {
+                break;
+            }
+        }
+        [x1, x2, x3]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::Distribution;
+    use super::UnitBall;
+
+    #[test]
+    fn value_stability() {
+        let mut rng = crate::test::rng(2);
+        let expected = [
+                [-0.42140140089381806, 0.4245276448803281, -0.7109276652167549],
+                [0.6683277779168173, 0.12753134283863998, 0.6843687153674809],
+                [-0.80397712218568, -0.0015797354643116712, 0.1588400395442835],
+            ];
+        let samples: [[f64; 3]; 3] = [
+                UnitBall.sample(&mut rng),
+                UnitBall.sample(&mut rng),
+                UnitBall.sample(&mut rng),
+            ];
+        assert_eq!(samples, expected);
+    }
+}

src/unit_disc.rs

@@ -0,0 +1,66 @@
+// Copyright 2019 Developers of the Rand project.
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
+// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
+// option. This file may not be copied, modified, or distributed
+// except according to those terms.
+
+use rand::Rng;
+use crate::{Distribution, Uniform, uniform::SampleUniform};
+use crate::utils::Float;
+
+/// Samples uniformly from the unit disc in two dimensions.
+///
+/// Implemented via rejection sampling.
+///
+///
+/// # Example
+///
+/// ```
+/// use rand_distr::{UnitDisc, Distribution};
+///
+/// let v: [f64; 2] = UnitDisc.sample(&mut rand::thread_rng());
+/// println!("{:?} is from the unit Disc.", v)
+/// ```
+#[derive(Clone, Copy, Debug)]
+pub struct UnitDisc;
+
+impl<N: Float + SampleUniform> Distribution<[N; 2]> for UnitDisc {
+    #[inline]
+    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [N; 2] {
+        let uniform = Uniform::new(N::from(-1.), N::from(1.));
+        let mut x1;
+        let mut x2;
+        loop {
+            x1 = uniform.sample(rng);
+            x2 = uniform.sample(rng);
+            if x1*x1 + x2*x2 <= N::from(1.) {
+                break;
+            }
+        }
+        [x1, x2]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::Distribution;
+    use super::UnitDisc;
+
+    #[test]
+    fn value_stability() {
+        let mut rng = crate::test::rng(2);
+        let expected = [
+            [-0.13921053103419823, -0.42140140089381806],
+            [0.4245276448803281, -0.7109276652167549],
+            [0.6683277779168173, 0.12753134283863998],
+        ];
+        let samples: [[f64; 2]; 3] = [
+                UnitDisc.sample(&mut rng),
+                UnitDisc.sample(&mut rng),
+                UnitDisc.sample(&mut rng),
+            ];
+        assert_eq!(samples, expected);
+    }
+}

src/unit_sphere.rs

@@ -1,4 +1,4 @@
-// Copyright 2018 Developers of the Rand project.
+// Copyright 2018-2019 Developers of the Rand project.
 //
 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
 // https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
@@ -18,19 +18,19 @@ use crate::utils::Float;
 /// # Example
 ///
 /// ```
-/// use rand_distr::{UnitSphereSurface, Distribution};
+/// use rand_distr::{UnitSphere, Distribution};
 ///
-/// let v: [f64; 3] = UnitSphereSurface.sample(&mut rand::thread_rng());
+/// let v: [f64; 3] = UnitSphere.sample(&mut rand::thread_rng());
 /// println!("{:?} is from the unit sphere surface.", v)
 /// ```
 ///
 /// [^1]: Marsaglia, George (1972). [*Choosing a Point from the Surface of a
 ///       Sphere.*](https://doi.org/10.1214/aoms/1177692644)
 ///       Ann. Math. Statist. 43, no. 2, 645--646.
 #[derive(Clone, Copy, Debug)]
-pub struct UnitSphereSurface;
+pub struct UnitSphere;
 
-impl<N: Float + SampleUniform> Distribution<[N; 3]> for UnitSphereSurface {
+impl<N: Float + SampleUniform> Distribution<[N; 3]> for UnitSphere {
     #[inline]
     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> [N; 3] {
         let uniform = Uniform::new(N::from(-1.), N::from(1.));
@@ -49,7 +49,7 @@ impl<N: Float + SampleUniform> Distribution<[N; 3]> for UnitSphereSurface {
 #[cfg(test)]
 mod tests {
     use crate::Distribution;
-    use super::UnitSphereSurface;
+    use super::UnitSphere;
 
     /// Assert that two numbers are almost equal to each other.
     ///
@@ -71,7 +71,7 @@ mod tests {
     fn norm() {
         let mut rng = crate::test::rng(1);
         for _ in 0..1000 {
-            let x: [f64; 3] = UnitSphereSurface.sample(&mut rng);
+            let x: [f64; 3] = UnitSphere.sample(&mut rng);
             assert_almost_eq!(x[0]*x[0] + x[1]*x[1] + x[2]*x[2], 1., 1e-15);
         }
     }
@@ -85,9 +85,9 @@ mod tests {
                 [0.9795722330927367, 0.18692349236651176, 0.07414747571708524],
             ];
         let samples: [[f64; 3]; 3] = [
-                UnitSphereSurface.sample(&mut rng),
-                UnitSphereSurface.sample(&mut rng),
-                UnitSphereSurface.sample(&mut rng),
+                UnitSphere.sample(&mut rng),
+                UnitSphere.sample(&mut rng),
+                UnitSphere.sample(&mut rng),
             ];
         assert_eq!(samples, expected);
     }

tests/uniformity.rs

@@ -6,16 +6,15 @@
 // option. This file may not be copied, modified, or distributed
 // except according to those terms.
 
-#![cfg(feature = "std")]
-
 #[macro_use]
 extern crate average;
 extern crate rand;
+extern crate rand_distr;
+extern crate core;
 
 use average::Histogram;
 use rand::distributions::Distribution;
 use rand::FromEntropy;
-use std as core;
 
 const N_BINS: usize = 100;
 const N_SAMPLES: u32 = 1_000_000;
@@ -28,10 +27,10 @@ fn unit_sphere() {
     const N_DIM: usize = 3;
     let h = Histogram100::with_const_width(-1., 1.);
     let mut histograms = [h.clone(), h.clone(), h];
-    let dist = rand::distributions::UnitSphereSurface::new();
+    let dist = rand_distr::UnitSphere;
     let mut rng = rand::rngs::SmallRng::from_entropy();
     for _ in 0..N_SAMPLES {
-        let v = dist.sample(&mut rng);
+        let v: [f64; 3] = dist.sample(&mut rng);
         for i in 0..N_DIM {
             histograms[i].add(v[i]).map_err(
                 |e| { println!("v: {}", v[i]); e }
@@ -52,10 +51,10 @@ fn unit_sphere() {
 fn unit_circle() {
     use std::f64::consts::PI;
     let mut h = Histogram100::with_const_width(-PI, PI);
-    let dist = rand::distributions::UnitCircle::new();
+    let dist = rand_distr::UnitCircle;
     let mut rng = rand::rngs::SmallRng::from_entropy();
     for _ in 0..N_SAMPLES {
-        let v = dist.sample(&mut rng);
+        let v: [f64; 2] = dist.sample(&mut rng);
         h.add(v[0].atan2(v[1])).unwrap();
     }
     let sum: u64 = h.bins().iter().sum();