lds_rs/lds.rs
1// #![feature(unboxed_closures)]
2//! Low-Discrepancy Sequence (LDS) Generator
3//!
4//! This code implements a set of low-discrepancy sequence generators, which are used to create sequences of numbers that are more evenly distributed than random numbers. These sequences are particularly useful in various fields such as computer graphics, numerical integration, and Monte Carlo simulations.
5//!
6//! The code defines several classes, each representing a different type of low-discrepancy sequence generator. The main types of sequences implemented are:
7//!
8//! 1. Van der Corput sequence
9//! 2. Halton sequence
10//! 3. Circle sequence
11//! 4. Sphere sequence
12//! 5. 3-Sphere Hopf sequence
13//! 6. N-dimensional Halton sequence
14//!
15//! Each generator takes specific inputs, usually in the form of base numbers or sequences of base numbers. These bases determine how the sequences are generated. The generators produce outputs in the form of floating-point numbers or lists of floating-point numbers, depending on the dimensionality of the sequence.
16//!
17//! The core algorithm used in most of these generators is the Van der Corput sequence. This sequence is created by expressing integers in a given base, reversing the digits, and placing them after a decimal point. For example, in base 2, the sequence would start: 1/2, 1/4, 3/4, 1/8, 5/8, and so on.
18//!
19//! The Halton sequence extends this concept to multiple dimensions by using a different base for each dimension. The Circle and Sphere sequences use trigonometric functions to map these low-discrepancy sequences onto circular or spherical surfaces.
20//!
21//! The code also includes utility functions and classes to support these generators. For instance, there's a list of prime numbers that can be used as bases for the sequences.
22//!
23//! Each generator class has methods to produce the next value in the sequence (pop()) and to reset the sequence to a specific starting point (reseed()). This allows for flexible use of the generators in various applications.
24//!
25//! The purpose of this code is to provide a toolkit for generating well-distributed sequences of numbers, which can be used in place of random numbers in many applications to achieve more uniform coverage of a given space or surface. This can lead to more efficient and accurate results in tasks like sampling, integration, and optimization.
26
27const TWO_PI: f64 = std::f64::consts::TAU;
28
29/// Van der Corput sequence
30///
31/// The `vdc` function is calculating the Van der Corput sequence value for a
32/// given index `k` and base `base`. It returns a `f64` value.
33///
34/// # Examples
35///
36/// ```rust
37/// use lds_rs::lds::vdc;
38///
39/// assert_eq!(vdc(11, 2), 0.8125);
40/// ```
41pub fn vdc(k: usize, base: usize) -> f64 {
42 let mut res = 0.0;
43 let mut denom = 1.0;
44 let mut k = k;
45 while k != 0 {
46 let remainder = k % base;
47 denom *= base as f64;
48 k /= base;
49 res += remainder as f64 / denom;
50 }
51 res
52}
53
54/// The `VdCorput` struct is a generator for the Van der Corput sequence, a low-discrepancy sequence
55/// commonly used in quasi-Monte Carlo methods.
56///
57/// Properties:
58///
59/// * `count`: The `count` property is used to keep track of the current iteration count of the Van der
60/// Corput sequence. It starts at 0 and increments by 1 each time the `pop()` method is called.
61/// * `base`: The `base` property represents the base of the Van der Corput sequence. It determines the
62/// number of digits used in each element of the sequence.
63///
64/// # Examples
65///
66/// ```rust
67/// use lds_rs::VdCorput;
68///
69/// let mut vgen = VdCorput::new(2);
70/// vgen.reseed(10);
71/// let result = vgen.pop();
72///
73/// assert_eq!(result, 0.8125);
74/// ```
75#[derive(Debug)]
76pub struct VdCorput {
77 count: usize,
78 base: usize,
79 rev_lst: [f64; 64],
80}
81
82impl VdCorput {
83 /// The `new` function creates a new [`VdCorput`] object with a given base for generating the Van der
84 /// Corput sequence.
85 ///
86 /// Arguments:
87 ///
88 /// * `base`: The `base` parameter is an integer value that is used to generate the Van der Corput
89 /// sequence. It determines the base of the sequence, which affects the distribution and pattern of the
90 /// generated numbers.
91 ///
92 /// Returns:
93 ///
94 /// The `new` function returns a `VdCorput` object.
95 pub fn new(base: usize) -> Self {
96 let mut rev_lst = [0.0; 64];
97 let mut reverse = 1.0;
98 for item in rev_lst.iter_mut() {
99 reverse /= base as f64;
100 *item = reverse;
101 }
102 VdCorput {
103 count: 0,
104 base,
105 rev_lst,
106 }
107 }
108
109 /// The `pop` function is a member function of the [`VdCorput`] class in Rust that increments the count
110 /// and calculates the next value in the Van der Corput sequence.
111 ///
112 /// Returns:
113 ///
114 /// The `pop` function returns a `f64` value, which is the next value in the Van der Corput sequence.
115 ///
116 /// # Examples
117 ///
118 /// ```rust
119 /// use lds_rs::lds::VdCorput;
120 ///
121 /// let mut vd_corput = VdCorput::new(2);
122 /// assert_eq!(vd_corput.pop(), 0.5);
123 /// ```
124 pub fn pop(&mut self) -> f64 {
125 self.count += 1; // ignore 0
126 let mut res = 0.0;
127 let mut k = self.count;
128 let mut i = 0;
129 while k != 0 {
130 let remainder = k % self.base;
131 k /= self.base;
132 res += remainder as f64 * self.rev_lst[i];
133 i += 1;
134 }
135 res
136 }
137
138 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
139 /// generator to a specific seed value. This allows the sequence generator to start generating the
140 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
141 /// seed.
142 pub fn reseed(&mut self, seed: usize) {
143 self.count = seed;
144 }
145}
146
147/// The [`Halton`] struct is a sequence generator that generates points in a 2-dimensional space using the
148/// Halton sequence.
149///
150/// Properties:
151///
152/// * `vdc0`: A variable of type [`VdCorput`] that represents the Van der Corput sequence generator for
153/// the first base. The Van der Corput sequence is a low-discrepancy sequence that is commonly used in
154/// quasi-Monte Carlo methods. It generates a sequence of numbers between 0 and
155/// * `vdc1`: The `vdc1` property is an instance of the [`VdCorput`] struct, which is responsible for
156/// generating the Van der Corput sequence with a base of 3. The Van der Corput sequence is another
157/// low-discrepancy sequence commonly used in quasi-Monte Carlo methods
158///
159/// # Examples
160///
161/// ```
162/// use lds_rs::Halton;
163///
164/// let mut hgen = Halton::new(&[2, 3]);
165/// hgen.reseed(10);
166/// let result = hgen.pop();
167/// assert_eq!(result[0], 0.8125);
168/// ```
169#[derive(Debug)]
170pub struct Halton {
171 vdc0: VdCorput,
172 vdc1: VdCorput,
173}
174
175impl Halton {
176 /// The `new` function creates a new [`Halton`] object with specified bases for generating the Halton
177 /// sequence.
178 ///
179 /// Arguments:
180 ///
181 /// * `base`: The `base` parameter is an array of two `usize` values. These values are used as the bases
182 /// for generating the Halton sequence. The first value in the array (`base[0]`) is used as the base for
183 /// generating the first component of the Halton sequence, and the second
184 ///
185 /// Returns:
186 ///
187 /// The `new` function returns an instance of the `Halton` struct.
188 pub fn new(base: &[usize]) -> Self {
189 Self {
190 vdc0: VdCorput::new(base[0]),
191 vdc1: VdCorput::new(base[1]),
192 }
193 }
194
195 /// Returns the pop of this [`Halton`].
196 ///
197 /// The `pop()` function is used to generate the next value in the sequence.
198 /// For example, in the [`VdCorput`] class, `pop()` increments the count and
199 /// calculates the Van der Corput sequence value for that count and base. In
200 /// the [`Halton`] class, `pop()` returns the next point in the Halton sequence
201 /// as a `[f64; 2]`. Similarly, in the `Circle` class, `pop()`
202 /// returns the next point on the unit circle as a `[f64; 2]`. In
203 /// the `Sphere` class, `pop()` returns the next point on the unit sphere as a
204 /// `[f64; 3]`. And in the `Sphere3Hopf` class, `pop()` returns
205 /// the next point on the 3-sphere using the Hopf fibration as a
206 /// `[f64; 4]`.
207 ///
208 /// Returns:
209 ///
210 /// An array of two f64 values is being returned.
211 ///
212 /// # Examples
213 ///
214 /// ```
215 /// use lds_rs::lds::Halton;
216 ///
217 /// let mut halton = Halton::new(&[2, 5]);
218 /// assert_eq!(halton.pop(), [0.5, 0.2]);
219 /// ```
220 pub fn pop(&mut self) -> [f64; 2] {
221 [self.vdc0.pop(), self.vdc1.pop()]
222 }
223
224 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
225 /// generator to a specific seed value. This allows the sequence generator to start generating the
226 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
227 /// seed.
228 #[allow(dead_code)]
229 pub fn reseed(&mut self, seed: usize) {
230 self.vdc0.reseed(seed);
231 self.vdc1.reseed(seed);
232 }
233}
234
235/// Circle sequence generator
236///
237/// The `Circle` struct is a generator for a circle sequence using the Van der Corput sequence.
238///
239/// Properties:
240///
241/// * `vdc`: A variable of type VdCorput, which is a sequence generator for Van der Corput sequence.
242///
243/// # Examples
244///
245/// ```
246/// use lds_rs::Circle;
247///
248/// let mut cgen = Circle::new(2);
249/// cgen.reseed(1);
250/// let result = cgen.pop();
251/// assert_eq!(result[1], 1.0);
252/// ```
253#[derive(Debug)]
254pub struct Circle {
255 vdc: VdCorput,
256}
257
258impl Circle {
259 /// Creates a new [`Circle`].
260 ///
261 /// The `new` function creates a new [`Circle`] object with a specified base value.
262 ///
263 /// Arguments:
264 ///
265 /// * `base`: The `base` parameter in the `new` function is the base value used to generate the Van
266 /// der Corput sequence. The Van der Corput sequence is a low-discrepancy sequence used in
267 /// quasi-Monte Carlo methods. It is generated by reversing the digits of the fractional part of the
268 ///
269 /// Returns:
270 ///
271 /// The `new` function is returning an instance of the `Circle` struct.
272 pub fn new(base: usize) -> Self {
273 Circle {
274 vdc: VdCorput::new(base),
275 }
276 }
277
278 /// Returns the pop of this [`Circle`].
279 ///
280 /// The `pop` function returns the coordinates of a point on a circle based on a random value.
281 ///
282 /// Returns:
283 ///
284 /// The `pop` function returns an array of two `f64` values, representing the sine and cosine of a
285 /// randomly generated angle.
286 ///
287 /// # Examples
288 ///
289 /// ```
290 /// use lds_rs::lds::Circle;
291 /// use approx_eq::assert_approx_eq;
292 ///
293 /// let mut circle = Circle::new(2);
294 /// let result = circle.pop();
295 /// assert_approx_eq!(result[1], 0.0);
296 /// assert_approx_eq!(result[0], -1.0);
297 /// ```
298 pub fn pop(&mut self) -> [f64; 2] {
299 // let two_pi = 2.0/// (-1.0 as f64).acos(); // ???
300 let theta = self.vdc.pop() * TWO_PI; // map to [0, 2*pi];
301 [theta.cos(), theta.sin()]
302 }
303
304 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
305 /// generator to a specific seed value. This allows the sequence generator to start generating the
306 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
307 /// seed.
308 #[allow(dead_code)]
309 pub fn reseed(&mut self, seed: usize) {
310 self.vdc.reseed(seed);
311 }
312}
313
314/// The [`Disk`] struct is a sequence generator that generates points in a 2-dimensional space using the
315/// Disk sequence.
316///
317/// Properties:
318///
319/// * `vdc0`: A variable of type [`VdCorput`] that represents the Van der Corput sequence generator for
320/// the first base. The Van der Corput sequence is a low-discrepancy sequence that is commonly used in
321/// quasi-Monte Carlo methods. It generates a sequence of numbers between 0 and
322/// * `vdc1`: The `vdc1` property is an instance of the [`VdCorput`] struct, which is responsible for
323/// generating the Van der Corput sequence with a base of 3. The Van der Corput sequence is another
324/// low-discrepancy sequence commonly used in quasi-Monte Carlo methods
325///
326/// # Examples
327///
328/// ```
329/// use lds_rs::Disk;
330/// use approx_eq::assert_approx_eq;
331///
332/// let mut dgen = Disk::new(&[2, 3]);
333/// dgen.reseed(0);
334/// let result = dgen.pop();
335/// assert_approx_eq!(result[0], -0.5773502691896257);
336/// ```
337#[derive(Debug)]
338pub struct Disk {
339 vdc0: VdCorput,
340 vdc1: VdCorput,
341}
342
343impl Disk {
344 /// The `new` function creates a new [`Disk`] object with specified bases for generating the Disk
345 /// sequence.
346 ///
347 /// Arguments:
348 ///
349 /// * `base`: The `base` parameter is an array of two `usize` values. These values are used as the bases
350 /// for generating the Disk sequence. The first value in the array (`base[0]`) is used as the base for
351 /// generating the first component of the Disk sequence, and the second
352 ///
353 /// Returns:
354 ///
355 /// The `new` function returns an instance of the `Disk` struct.
356 pub fn new(base: &[usize]) -> Self {
357 Self {
358 vdc0: VdCorput::new(base[0]),
359 vdc1: VdCorput::new(base[1]),
360 }
361 }
362
363 /// Returns the pop of this [`Disk`].
364 ///
365 /// The `pop()` function is used to generate the next value in the sequence.
366 /// For example, in the [`VdCorput`] class, `pop()` increments the count and
367 /// calculates the Van der Corput sequence value for that count and base. In
368 /// the [`Disk`] class, `pop()` returns the next point in the Disk sequence
369 /// as a `[f64; 2]`. Similarly, in the `Circle` class, `pop()`
370 /// returns the next point on the unit circle as a `[f64; 2]`. In
371 /// the `Sphere` class, `pop()` returns the next point on the unit sphere as a
372 /// `[f64; 3]`. And in the `Sphere3Hopf` class, `pop()` returns
373 /// the next point on the 3-sphere using the Hopf fibration as a
374 /// `[f64; 4]`.
375 ///
376 /// Returns:
377 ///
378 /// An array of two f64 values is being returned.
379 ///
380 /// # Examples
381 ///
382 /// ```
383 /// use lds_rs::lds::Disk;
384 /// use approx_eq::assert_approx_eq;
385 ///
386 /// let mut dgen = Disk::new(&[2, 3]);
387 /// let result = dgen.pop();
388 /// assert_approx_eq!(result[0], -0.5773502691896257);
389 /// ```
390 pub fn pop(&mut self) -> [f64; 2] {
391 let theta = self.vdc0.pop() * TWO_PI; // map to [0, 2*pi];
392 let radius = self.vdc1.pop().sqrt();
393 [radius * theta.cos(), radius * theta.sin()]
394 }
395
396 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
397 /// generator to a specific seed value. This allows the sequence generator to start generating the
398 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
399 /// seed.
400 #[allow(dead_code)]
401 pub fn reseed(&mut self, seed: usize) {
402 self.vdc0.reseed(seed);
403 self.vdc1.reseed(seed);
404 }
405}
406
407/// Sphere sequence generator
408///
409/// The `Sphere` struct is a generator for a sequence of points on a sphere.
410///
411/// Properties:
412///
413/// * `vdc`: The `vdc` property is an instance of the [`VdCorput`] struct. It is used to generate a Van
414/// der Corput sequence, which is a low-discrepancy sequence used for sampling points in a unit
415/// interval.
416/// * `cirgen`: The `cirgen` property is an instance of the [`Circle`] struct. It is responsible for
417/// generating points on a circle.
418///
419/// # Examples
420///
421/// ```
422/// use lds_rs::Sphere;
423///
424/// let mut sgen = Sphere::new(&[2, 3]);
425/// sgen.reseed(1);
426/// let result = sgen.pop();
427/// assert_eq!(result[2], -0.5);
428/// ```
429#[derive(Debug)]
430pub struct Sphere {
431 vdc: VdCorput,
432 cirgen: Circle,
433}
434
435impl Sphere {
436 /// Creates a new [`Sphere`].
437 ///
438 /// The function `new` creates a new [`Sphere`] object with a given base.
439 ///
440 /// Arguments:
441 ///
442 /// * `base`: The `base` parameter is an array of `usize` values. It is used to initialize the `Sphere`
443 /// struct. The first element of the `base` array is used to create a new `VdCorput` struct, and the
444 /// second element is used to create a new `Circle
445 ///
446 /// Returns:
447 ///
448 /// The `new` function is returning an instance of the `Sphere` struct.
449 pub fn new(base: &[usize]) -> Self {
450 Sphere {
451 vdc: VdCorput::new(base[0]),
452 cirgen: Circle::new(base[1]),
453 }
454 }
455
456 /// Returns the pop of this [`Sphere`].
457 ///
458 /// The `pop` function returns a random point on a sphere using the VDC and cirgen generators.
459 ///
460 /// Returns:
461 ///
462 /// an array of three `f64` values, representing the coordinates of a point on a sphere. The first
463 /// two values (`sinphi * c` and `sinphi * s`) represent the x and y coordinates, while the third
464 /// value (`cosphi`) represents the z coordinate.
465 ///
466 /// # Examples
467 ///
468 /// ```
469 /// use lds_rs::lds::Sphere;
470 /// use approx_eq::assert_approx_eq;
471 ///
472 /// let mut sphere = Sphere::new(&[2, 3]);
473 /// let result = sphere.pop();
474 /// assert_approx_eq!(result[0], -0.5);
475 /// assert_approx_eq!(result[1], 0.8660254037844387);
476 /// assert_approx_eq!(result[2], 0.0);
477 /// ```
478 pub fn pop(&mut self) -> [f64; 3] {
479 let cosphi = 2.0 * self.vdc.pop() - 1.0; // map to [-1, 1];
480 let sinphi = (1.0 - cosphi * cosphi).sqrt();
481 let [c, s] = self.cirgen.pop();
482 [sinphi * c, sinphi * s, cosphi]
483 }
484
485 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
486 /// generator to a specific seed value. This allows the sequence generator to start generating the
487 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
488 /// seed.
489 #[allow(dead_code)]
490 pub fn reseed(&mut self, seed: usize) {
491 self.cirgen.reseed(seed);
492 self.vdc.reseed(seed);
493 }
494}
495
496/// The `Sphere3Hopf` struct is a sequence generator for the S(3) sequence using Hopf coordinates.
497///
498/// Properties:
499///
500/// * `vdc0`: An instance of the VdCorput sequence generator used for the first coordinate of the Hopf
501/// coordinates.
502/// * `vdc1`: The `vdc1` property is an instance of the [`VdCorput`] struct, which is used to generate a
503/// Van der Corput sequence. This sequence is a low-discrepancy sequence that is commonly used in
504/// numerical methods for generating random numbers. In this case, it is
505/// * `vdc2`: The `vdc2` property is an instance of the [`VdCorput`] struct, which is used to generate a
506/// Van der Corput sequence. This sequence is a low-discrepancy sequence that is commonly used in
507/// numerical methods for generating random numbers. In the context of the `
508///
509/// The `Sphere3Hopf` class is a sequence generator that generates points on a
510/// 3-sphere using the Hopf fibration. It uses three instances of the `VdCorput`
511/// class to generate the sequence values and maps them to points on the
512/// 3-sphere. The `pop()` method returns the next point on the 3-sphere as a
513/// `[f64; 4]`, where the first three elements represent the x, y,
514/// and z coordinates of the point, and the fourth element represents the w
515/// coordinate. The `reseed()` method is used to reset the state of the sequence
516/// generator to a specific seed value.
517///
518/// # Examples
519///
520/// ```
521/// use lds_rs::Sphere3Hopf;
522/// use approx_eq::assert_approx_eq;
523///
524/// let mut sgen = Sphere3Hopf::new(&[2, 3, 5]);
525/// sgen.reseed(0);
526/// let result = sgen.pop();
527/// assert_approx_eq!(result[2], 0.4472135954999573);
528/// ```
529#[derive(Debug)]
530pub struct Sphere3Hopf {
531 vdc0: VdCorput,
532 vdc1: VdCorput,
533 vdc2: VdCorput,
534}
535
536impl Sphere3Hopf {
537 /// Creates a new [`Sphere3Hopf`].
538 ///
539 /// The `new` function creates a new instance of the [`Sphere3Hopf`] struct with three `VdCorput`
540 /// instances initialized with the values from the `base` slice.
541 ///
542 /// Arguments:
543 ///
544 /// * `base`: The `base` parameter is an array of three `usize` values. These values are used to
545 /// initialize three instances of the `VdCorput` struct, which is a type of quasi-random number
546 /// generator. Each `VdCorput` instance is initialized with a different base value from the
547 ///
548 /// Returns:
549 ///
550 /// The `new` function is returning an instance of the `Sphere3Hopf` struct.
551 pub fn new(base: &[usize]) -> Self {
552 Sphere3Hopf {
553 vdc0: VdCorput::new(base[0]),
554 vdc1: VdCorput::new(base[1]),
555 vdc2: VdCorput::new(base[2]),
556 }
557 }
558
559 /// The `pop` function returns a four-element array representing the coordinates of a point on a
560 /// sphere in 3D space.
561 ///
562 /// Returns:
563 ///
564 /// The function `pop` returns an array of four `f64` values.
565 /// Returns the pop of this [`Sphere3Hopf`].
566 ///
567 /// The `pop()` function is used to generate the next value in the sequence.
568 /// For example, in the [`VdCorput`] class, `pop()` increments the count and
569 /// calculates the Van der Corput sequence value for that count and base. In
570 /// the [`Disk`] class, `pop()` returns the next point in the Disk sequence
571 /// as a `[f64; 2]`. Similarly, in the [`Circle`] class, `pop()`
572 /// returns the next point on the unit circle as a `[f64; 2]`. In
573 /// the [`Sphere`] class, `pop()` returns the next point on the unit sphere as a
574 /// `[f64; 3]`. And in the [`Sphere3Hopf`] class, `pop()` returns
575 /// the next point on the 3-sphere using the Hopf fibration as a `[f64; 4]`.
576 ///
577 /// # Examples
578 ///
579 /// ```
580 /// use lds_rs::Sphere3Hopf;
581 /// use approx_eq::assert_approx_eq;
582 ///
583 /// let mut sgen = Sphere3Hopf::new(&[2, 3, 5]);
584 /// let result = sgen.pop();
585 /// assert_approx_eq!(result[2], 0.4472135954999573);
586 pub fn pop(&mut self) -> [f64; 4] {
587 let phi = self.vdc0.pop() * TWO_PI; // map to [0, 2*pi];
588 let psy = self.vdc1.pop() * TWO_PI; // map to [0, 2*pi];
589 let vdc = self.vdc2.pop();
590 let cos_eta = vdc.sqrt();
591 let sin_eta = (1.0 - vdc).sqrt();
592 [
593 cos_eta * psy.cos(),
594 cos_eta * psy.sin(),
595 sin_eta * (phi + psy).cos(),
596 sin_eta * (phi + psy).sin(),
597 ]
598 }
599
600 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
601 /// generator to a specific seed value. This allows the sequence generator to start generating the
602 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
603 /// seed.
604 #[allow(dead_code)]
605 pub fn reseed(&mut self, seed: usize) {
606 self.vdc0.reseed(seed);
607 self.vdc1.reseed(seed);
608 self.vdc2.reseed(seed);
609 }
610}
611
612// use crate::lds::VdCorput;
613
614/// The HaltonN struct is a generator for the Halton(n) sequence.
615///
616/// Properties:
617///
618/// * `vdcs`: A vector of VdCorput objects.
619#[derive(Debug)]
620pub struct HaltonN {
621 vdcs: Vec<VdCorput>,
622}
623
624/// Halton(n) sequence generator
625///
626/// The [`HaltonN`] class is a sequence generator that generates points in a
627/// N-dimensional space using the Halton sequence. The Halton sequence is a
628/// low-discrepancy sequence that is commonly used in quasi-Monte Carlo methods.
629/// It is generated by iterating over two different bases and calculating the
630/// fractional parts of the numbers in those bases. The [`HaltonN`] class keeps
631/// track of the current count and bases, and provides a `pop()` method that
632/// returns the next point in the sequence as a `std::vector<double>`.
633///
634/// # Examples
635///
636/// ```
637/// use lds_rs::HaltonN;
638/// use approx_eq::assert_approx_eq;
639///
640/// let mut hgen = HaltonN::new(&[2, 3, 5, 7, 11]);
641/// hgen.reseed(10);
642/// for _i in 0..10 {
643/// println!("{:?}", hgen.pop_vec());
644/// }
645/// let res = hgen.pop_vec();
646///
647/// assert_approx_eq!(res[0], 0.65625);
648/// ```
649impl HaltonN {
650 /// Creates a new [`HaltonN`].
651 ///
652 /// The `new` function creates a new `HaltonN` struct with a specified number of `VdCorput`
653 /// instances.
654 ///
655 /// Arguments:
656 ///
657 /// * `base`: The `base` parameter is a slice of `usize` values. It represents the base values for
658 /// each dimension of the Halton sequence. Each dimension of the Halton sequence uses a different
659 /// base value to generate the sequence.
660 ///
661 /// Returns:
662 ///
663 /// The `new` function returns a new instance of the `HaltonN` struct.
664 pub fn new(base: &[usize]) -> Self {
665 HaltonN {
666 vdcs: base.iter().map(|b| VdCorput::new(*b)).collect(),
667 // vdcs: (0..n).map(|i| VdCorput::new(base[i])).collect(),
668 }
669 }
670
671 /// Returns the pop vec of this [`HaltonN`].
672 ///
673 /// The `pop()` function is used to generate the next value in the sequence.
674 /// For example, in the [`VdCorput`] class, `pop()` increments the count and
675 /// calculates the Van der Corput sequence value for that count and base. In
676 /// the [`Halton`] class, `pop()` returns the next point in the Halton sequence
677 /// as a `[f64; 2]`. Similarly, in the [`Circle`] class, `pop()`
678 /// returns the next point on the unit circle as a `[f64; 2]`. In
679 /// the [`Sphere`] class, `pop()` returns the next point on the unit sphere as a
680 /// `[f64; 3]`. And in the [`Sphere3Hopf`] class, `pop()` returns
681 /// the next point on the 3-sphere using the Hopf fibration as a `[f64; 4]`.
682 pub fn pop_vec(&mut self) -> Vec<f64> {
683 self.vdcs.iter_mut().map(|vdc| vdc.pop()).collect()
684 }
685
686 /// The below code is a Rust function called `reseed` that is used to reset the state of a sequence
687 /// generator to a specific seed value. This allows the sequence generator to start generating the
688 /// sequence from the beginning or from a specific point in the sequence, depending on the value of the
689 /// seed.
690 pub fn reseed(&mut self, seed: usize) {
691 for vdc in self.vdcs.iter_mut() {
692 vdc.reseed(seed);
693 }
694 }
695}
696
697// First 1000 prime numbers;
698#[allow(dead_code)]
699pub const PRIME_TABLE: [usize; 1000] = [
700 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
701 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193,
702 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
703 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421,
704 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547,
705 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
706 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
707 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929,
708 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039,
709 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153,
710 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
711 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409,
712 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
713 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613,
714 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
715 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873,
716 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999,
717 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113,
718 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251,
719 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371,
720 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477,
721 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647,
722 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731,
723 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857,
724 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001,
725 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163,
726 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299,
727 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407,
728 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539,
729 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659,
730 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793,
731 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919,
732 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051,
733 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201,
734 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327,
735 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463,
736 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603,
737 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733,
738 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903,
739 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009,
740 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153,
741 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303,
742 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441,
743 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569,
744 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701,
745 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843,
746 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987,
747 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131,
748 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269,
749 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373,
750 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553,
751 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691,
752 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829,
753 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967,
754 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109,
755 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247,
756 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451,
757 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559,
758 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687,
759 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841,
760 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
761];
762
763#[cfg(test)]
764mod tests {
765 use super::*;
766 use approx_eq::assert_approx_eq;
767
768 #[test]
769 fn test_vdc() {
770 assert_approx_eq!(vdc(11, 2), 0.8125);
771 }
772
773 #[test]
774 fn test_vdcorput() {
775 let mut vgen = VdCorput::new(2);
776 vgen.reseed(0);
777 assert_approx_eq!(vgen.pop(), 0.5);
778 assert_approx_eq!(vgen.pop(), 0.25);
779 assert_approx_eq!(vgen.pop(), 0.75);
780
781 let mut vgen = VdCorput::new(3);
782 vgen.reseed(0);
783 assert_approx_eq!(vgen.pop(), 1.0/3.0);
784 assert_approx_eq!(vgen.pop(), 2.0/3.0);
785 assert_approx_eq!(vgen.pop(), 1.0/9.0);
786 }
787
788 #[test]
789 fn test_halton() {
790 let mut hgen = Halton::new(&[2, 3]);
791 hgen.reseed(0);
792 let res = hgen.pop();
793 assert_approx_eq!(res[0], 0.5);
794 assert_approx_eq!(res[1], 1.0/3.0);
795 let res = hgen.pop();
796 assert_approx_eq!(res[0], 0.25);
797 assert_approx_eq!(res[1], 2.0/3.0);
798 }
799
800 #[test]
801 fn test_circle() {
802 let mut cgen = Circle::new(2);
803 cgen.reseed(0);
804 let res = cgen.pop();
805 assert_approx_eq!(res[0], -1.0);
806 assert_approx_eq!(res[1], 0.0);
807 let res = cgen.pop();
808 assert_approx_eq!(res[0], 0.0);
809 assert_approx_eq!(res[1], 1.0);
810 }
811
812 #[test]
813 fn test_disk() {
814 let mut dgen = Disk::new(&[2, 3]);
815 dgen.reseed(0);
816 let res = dgen.pop();
817 assert_approx_eq!(res[0], -0.5773502691896258);
818 assert_approx_eq!(res[1], 0.0);
819 }
820
821 #[test]
822 fn test_sphere() {
823 let mut sgen = Sphere::new(&[2, 3]);
824 sgen.reseed(0);
825 let res = sgen.pop();
826 assert_approx_eq!(res[0], -0.5);
827 assert_approx_eq!(res[1], 0.8660254037844387);
828 assert_approx_eq!(res[2], 0.0);
829 let res = sgen.pop();
830 assert_approx_eq!(res[0], -0.4330127018922193);
831 assert_approx_eq!(res[1], -0.75);
832 assert_approx_eq!(res[2], -0.5);
833 }
834
835 #[test]
836 fn test_sphere3hopf() {
837 let mut sgen = Sphere3Hopf::new(&[2, 3, 5]);
838 sgen.reseed(0);
839 let res = sgen.pop();
840 assert_approx_eq!(res[0], -0.22360679774997885);
841 assert_approx_eq!(res[1], 0.3872983346207417);
842 assert_approx_eq!(res[2], 0.44721359549995726);
843 assert_approx_eq!(res[3], -0.7745966692414837);
844 }
845
846 #[test]
847 fn test_halton_n() {
848 let mut hgen = HaltonN::new(&[2, 3, 5]);
849 hgen.reseed(0);
850 let res = hgen.pop_vec();
851 assert_approx_eq!(res[0], 0.5);
852 assert_approx_eq!(res[1], 1.0/3.0);
853 assert_approx_eq!(res[2], 0.2);
854 let res = hgen.pop_vec();
855 assert_approx_eq!(res[0], 0.25);
856 assert_approx_eq!(res[1], 2.0/3.0);
857 assert_approx_eq!(res[2], 0.4);
858 }
859}