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}