egobox_doe/
lhs.rs

1use crate::SamplingMethod;
2use crate::utils::{cdist, pdist};
3use linfa::Float;
4use ndarray::{Array, Array2, ArrayBase, Axis, Data, Ix2, ShapeBuilder, s};
5use ndarray_rand::{
6    RandomExt, rand::Rng, rand::SeedableRng, rand::seq::SliceRandom, rand_distr::Uniform,
7};
8use ndarray_stats::QuantileExt;
9use rand_xoshiro::Xoshiro256Plus;
10use std::cmp;
11use std::sync::{Arc, RwLock};
12
13#[cfg(feature = "serializable")]
14use serde::{Deserialize, Serialize};
15
16/// Kinds of Latin Hypercube Design
17#[derive(Clone, Debug, Default, Copy)]
18#[cfg_attr(feature = "serializable", derive(Serialize, Deserialize))]
19pub enum LhsKind {
20    /// sample is choosen randomly within its latin hypercube intervals
21    Classic,
22    /// sample is the middle of its latin hypercube intervals
23    Centered,
24    /// distance between points is maximized
25    Maximin,
26    /// sample is the middle of its latin hypercube intervals and distance between points is maximized
27    CenteredMaximin,
28    /// samples locations is optimized using the Enhanced Stochastic Evolutionary algorithm (ESE)
29    /// See Jin, R. and Chen, W. and Sudjianto, A. (2005), “An efficient algorithm for constructing
30    /// optimal design of computer experiments.” Journal of Statistical Planning and Inference, 134:268-287.
31    #[default]
32    Optimized,
33}
34
35type RngRef<R> = Arc<RwLock<R>>;
36
37/// The LHS design is built as follows: each dimension space is divided into ns sections
38/// where ns is the number of sampling points, and one point in selected in each section.
39/// The selection method gives different kind of LHS (see [LhsKind])
40#[derive(Clone, Debug)]
41#[cfg_attr(feature = "serializable", derive(Serialize, Deserialize))]
42pub struct Lhs<F: Float, R: Rng> {
43    /// Sampling space definition as a (nx, 2) matrix
44    /// The ith row is the [lower_bound, upper_bound] of xi, the ith component of x
45    xlimits: Array2<F>,
46    /// The requested kind of LHS
47    kind: LhsKind,
48    /// Random generator used for reproducibility (not used in case of Centered LHS)
49    rng: RngRef<R>,
50}
51
52/// LHS with default random generator
53impl<F: Float> Lhs<F, Xoshiro256Plus> {
54    /// Constructor given a design space given a (nx, 2) matrix \[\[lower bound, upper bound\], ...\]
55    ///
56    /// ```
57    /// use egobox_doe::Lhs;
58    /// use ndarray::arr2;
59    ///
60    /// let doe = Lhs::new(&arr2(&[[0.0, 1.0], [5.0, 10.0]]));
61    /// ```
62    pub fn new(xlimits: &ArrayBase<impl Data<Elem = F>, Ix2>) -> Self {
63        Self::new_with_rng(xlimits, Xoshiro256Plus::from_entropy())
64    }
65}
66
67impl<F: Float, R: Rng> SamplingMethod<F> for Lhs<F, R> {
68    fn sampling_space(&self) -> &Array2<F> {
69        &self.xlimits
70    }
71
72    fn normalized_sample(&self, ns: usize) -> Array2<F> {
73        match &self.kind {
74            LhsKind::Classic => self._classic_lhs(ns),
75            LhsKind::Centered => self._centered_lhs(ns),
76            LhsKind::Maximin => self._maximin_lhs(ns, false, 5),
77            LhsKind::CenteredMaximin => self._maximin_lhs(ns, true, 5),
78            LhsKind::Optimized => {
79                let doe = self._classic_lhs(ns);
80                let nx = self.xlimits.nrows();
81                let outer_loop = cmp::min((1.5 * nx as f64) as usize, 30);
82                let inner_loop = cmp::min(20 * nx, 100);
83                self._maximin_ese(&doe, outer_loop, inner_loop)
84            }
85        }
86    }
87}
88
89impl<F: Float, R: Rng> Lhs<F, R> {
90    /// Constructor with given design space and random generator.
91    /// * `xlimits`: (nx, 2) matrix where nx is the dimension of the samples and the ith row
92    ///   is the definition interval of the ith component of x.
93    /// * `rng`: random generator used for [LhsKind::Classic] and [LhsKind::Optimized] LHS
94    pub fn new_with_rng(xlimits: &ArrayBase<impl Data<Elem = F>, Ix2>, rng: R) -> Self {
95        if xlimits.ncols() != 2 {
96            panic!("xlimits must have 2 columns (lower, upper)");
97        }
98        Lhs {
99            xlimits: xlimits.to_owned(),
100            kind: LhsKind::default(),
101            rng: Arc::new(RwLock::new(rng)),
102        }
103    }
104
105    /// Sets the kind of LHS
106    pub fn kind(mut self, kind: LhsKind) -> Self {
107        self.kind = kind;
108        self
109    }
110
111    /// Sets the random generator
112    pub fn with_rng<R2: Rng>(self, rng: R2) -> Lhs<F, R2> {
113        Lhs {
114            xlimits: self.xlimits,
115            kind: self.kind,
116            rng: Arc::new(RwLock::new(rng)),
117        }
118    }
119
120    fn _maximin_ese(&self, lhs: &Array2<F>, outer_loop: usize, inner_loop: usize) -> Array2<F> {
121        // hard-coded params
122        let j_range = 20;
123        let p = F::cast(10.);
124        let t0 = F::cast(0.005) * self._phip(lhs, p);
125        let tol = F::cast(1e-3);
126
127        let mut t = t0;
128        let mut lhs_own = lhs.to_owned();
129        let mut lhs_best = lhs.to_owned();
130        let nx = lhs.ncols();
131        let mut phip = self._phip(&lhs_best, p);
132        let mut phip_best = phip;
133
134        for _ in 0..outer_loop {
135            let mut n_acpt = 0.;
136            let mut n_imp = 0.;
137
138            for i in 0..inner_loop {
139                let modulo = (i + 1) % nx;
140                let mut l_x: Vec<Array2<F>> = Vec::with_capacity(j_range);
141                let mut l_phip: Vec<F> = Vec::with_capacity(j_range);
142
143                // Build j different plans with a single swap procedure
144                // See description of phip_swap procedure
145                let mut rng = self.rng.write().unwrap();
146                for j in 0..j_range {
147                    l_x.push(lhs_own.to_owned());
148                    let php = self._phip_swap(&mut l_x[j], modulo, phip, p, &mut *rng);
149                    l_phip.push(php);
150                }
151                let lphip = Array::from_shape_vec(j_range, l_phip).unwrap();
152                let k = lphip.argmin().unwrap();
153                let phip_try = lphip[k];
154                // Threshold of acceptance
155                if phip_try - phip <= t * F::cast(rng.r#gen::<f64>()) {
156                    phip = phip_try;
157                    n_acpt += 1.;
158                    lhs_own = l_x.swap_remove(k);
159
160                    // best plan retained
161                    if phip < phip_best {
162                        lhs_best = lhs_own.to_owned();
163                        phip_best = phip;
164                        n_imp += 1.;
165                    }
166                }
167            }
168            let p_accpt = n_acpt / (inner_loop as f64); // probability of acceptance
169            let p_imp = n_imp / (inner_loop as f64); // probability of improvement
170
171            if phip - phip_best > tol {
172                if p_accpt >= 0.1 && p_imp < p_accpt {
173                    t *= F::cast(0.8)
174                } else if p_accpt >= 0.1 && (p_imp - p_accpt).abs() < f64::EPSILON {
175                } else {
176                    t /= F::cast(0.8)
177                }
178            } else if p_accpt <= 0.1 {
179                t /= F::cast(0.7)
180            } else {
181                t *= F::cast(0.9)
182            }
183        }
184        lhs_best
185    }
186
187    fn _phip(&self, lhs: &ArrayBase<impl Data<Elem = F>, Ix2>, p: F) -> F {
188        F::powf(pdist(lhs).mapv(|v| F::powf(v, -p)).sum(), F::one() / p)
189    }
190
191    fn _phip_swap(&self, x: &mut Array2<F>, k: usize, phip: F, p: F, rng: &mut R) -> F {
192        // Choose two random rows
193        let i1 = rng.gen_range(0..x.nrows());
194        let mut i2 = rng.gen_range(0..x.nrows());
195        while i2 == i1 {
196            i2 = rng.gen_range(0..x.nrows());
197        }
198        // Compute new phip
199        let mut x_rest = Array2::zeros((x.nrows() - 2, x.ncols()));
200        let mut row_i = 0;
201        for (i, row) in x.axis_iter(Axis(0)).enumerate() {
202            if i != i1 && i != i2 {
203                x_rest.row_mut(row_i).assign(&row);
204                row_i += 1;
205            }
206        }
207
208        let mut dist1 = cdist(&x.slice(s![i1..i1 + 1, ..]), &x_rest);
209        let mut dist2 = cdist(&x.slice(s![i2..i2 + 1, ..]), &x_rest);
210
211        let m1 = (x_rest.column(k).to_owned() - x[[i1, k]]).map(|v| *v * *v);
212        let m2 = (x_rest.column(k).to_owned() - x[[i2, k]]).map(|v| *v * *v);
213
214        let two = F::cast(2.);
215        let mut d1 = dist1.mapv(|v| v * v) - &m1 + &m2;
216        d1.mapv_inplace(|v| F::powf(v, -p / two));
217        let mut d2 = dist2.mapv(|v| v * v) + &m1 - &m2;
218        d2.mapv_inplace(|v| F::powf(v, -p / two));
219
220        dist1.mapv_inplace(|v| F::powf(v, -p));
221        dist2.mapv_inplace(|v| F::powf(v, -p));
222        let mut res = (d1 - dist1).sum() + (d2 - dist2).sum();
223        res = F::powf(F::powf(phip, p) + res, F::one() / p);
224
225        // swap points
226        x.swap([i1, k], [i2, k]);
227        res
228    }
229
230    fn _classic_lhs(&self, ns: usize) -> Array2<F> {
231        let nx = self.xlimits.nrows();
232        let cut = Array::linspace(0., 1., ns + 1);
233
234        let mut rng = self.rng.write().unwrap();
235        let rnd = Array::random_using((ns, nx).f(), Uniform::new(0., 1.), &mut *rng);
236        let a = cut.slice(s![..ns]).to_owned();
237        let b = cut.slice(s![1..(ns + 1)]);
238        let c = &b - &a;
239        let mut rdpoints = Array::zeros((ns, nx).f());
240        for j in 0..nx {
241            let d = rnd.column(j).to_owned() * &c + &a;
242            rdpoints.column_mut(j).assign(&d)
243        }
244        let mut lhs = Array::zeros((ns, nx).f());
245        for j in 0..nx {
246            let mut colj = rdpoints.column_mut(j);
247            colj.as_slice_mut().unwrap().shuffle(&mut *rng);
248            lhs.column_mut(j).assign(&colj);
249        }
250        lhs.mapv_into_any(F::cast)
251    }
252
253    fn _centered_lhs(&self, ns: usize) -> Array2<F> {
254        let nx = self.xlimits.nrows();
255        let cut = Array::linspace(0., 1., ns + 1);
256
257        let a = cut.slice(s![..ns]).to_owned();
258        let b = cut.slice(s![1..(ns + 1)]);
259        let mut c = (a + b) / 2.;
260        let mut lhs = Array::zeros((ns, nx).f());
261
262        let mut rng = self.rng.write().unwrap();
263        for j in 0..nx {
264            c.as_slice_mut().unwrap().shuffle(&mut *rng);
265            lhs.column_mut(j).assign(&c);
266        }
267        lhs.mapv_into_any(F::cast)
268    }
269
270    fn _maximin_lhs(&self, ns: usize, centered: bool, max_iters: usize) -> Array2<F> {
271        let mut lhs = if centered {
272            self._centered_lhs(ns)
273        } else {
274            self._classic_lhs(ns)
275        };
276        let mut max_dist = *pdist(&lhs).min().unwrap();
277        let mut lhs_maximin = lhs;
278        for _ in 0..max_iters - 1 {
279            if centered {
280                lhs = self._centered_lhs(ns);
281            } else {
282                lhs = self._classic_lhs(ns);
283            }
284            let d_min = *pdist(&lhs).min().unwrap();
285            if max_dist < d_min {
286                max_dist = d_min;
287                std::mem::swap(&mut lhs_maximin, &mut lhs)
288            }
289        }
290        lhs_maximin
291    }
292}
293
294#[cfg(test)]
295mod tests {
296    use super::*;
297    use approx::{assert_abs_diff_eq, assert_abs_diff_ne};
298    use ndarray::{arr2, array};
299    use std::time::Instant;
300
301    #[test]
302    fn test_lhs() {
303        let xlimits = arr2(&[[5., 10.], [0., 1.]]);
304        let expected = array![
305            [9.000042958859238, 0.2175219214807449],
306            [5.085755595295461, 0.7725590934255249],
307            [7.062569781563214, 0.44540674774531397],
308            [8.306461322653673, 0.9046507902710129],
309            [6.310411395727105, 0.0606130622609971]
310        ];
311        let actual = Lhs::new(&xlimits)
312            .with_rng(Xoshiro256Plus::seed_from_u64(42))
313            .sample(5);
314        assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
315    }
316
317    #[test]
318    fn test_lhs_speed() {
319        let start = Instant::now();
320        let xlimits = arr2(&[[0., 1.], [0., 1.]]);
321        let n = 10;
322        let _actual = Lhs::new(&xlimits).sample(n);
323        let duration = start.elapsed();
324        println!("Time elapsed in optimized LHS is: {duration:?}");
325    }
326
327    #[test]
328    fn test_classic_lhs() {
329        let xlimits = arr2(&[[5., 10.], [0., 1.]]);
330        let expected = array![
331            [9.000042958859238, 0.44540674774531397],
332            [5.085755595295461, 0.7725590934255249],
333            [7.062569781563214, 0.2175219214807449],
334            [8.306461322653673, 0.9046507902710129],
335            [6.310411395727105, 0.0606130622609971]
336        ];
337        let actual = Lhs::new(&xlimits)
338            .with_rng(Xoshiro256Plus::seed_from_u64(42))
339            .kind(LhsKind::Classic)
340            .sample(5);
341        assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
342    }
343
344    #[test]
345    fn test_centered_lhs() {
346        let xlimits = arr2(&[[5., 10.], [0., 1.]]);
347        let expected = array![[7.5, 0.9], [8.5, 0.1], [5.5, 0.7], [6.5, 0.3], [9.5, 0.5]];
348        let actual = Lhs::new(&xlimits)
349            .with_rng(Xoshiro256Plus::seed_from_u64(0))
350            .kind(LhsKind::Centered)
351            .sample(5);
352        assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
353    }
354
355    #[test]
356    fn test_centered_maximin_lhs() {
357        let xlimits = arr2(&[[5., 10.], [0., 1.]]);
358        let expected = array![[7.5, 0.9], [8.5, 0.1], [5.5, 0.7], [6.5, 0.3], [9.5, 0.5]];
359        let actual = Lhs::new(&xlimits)
360            .with_rng(Xoshiro256Plus::seed_from_u64(0))
361            .kind(LhsKind::CenteredMaximin)
362            .sample(5);
363        assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
364    }
365
366    #[test]
367    fn test_phip_swap() {
368        let xlimits = arr2(&[[0., 1.], [0., 1.]]);
369        let k = 1;
370        let phip = 7.290525742903316;
371        let mut p0 = array![
372            [0.45, 0.75],
373            [0.75, 0.95],
374            [0.05, 0.45],
375            [0.55, 0.15000000000000002],
376            [0.35000000000000003, 0.25],
377            [0.95, 0.8500000000000001],
378            [0.15000000000000002, 0.55],
379            [0.25, 0.05],
380            [0.8500000000000001, 0.35000000000000003],
381            [0.6500000000000001, 0.6500000000000001]
382        ];
383        let p = 10.;
384        let mut rng = Xoshiro256Plus::seed_from_u64(42);
385        let _res = Lhs::new(&xlimits)._phip_swap(&mut p0, k, phip, p, &mut rng);
386    }
387
388    #[test]
389    fn test_no_duplicate() {
390        let xlimits = arr2(&[[5., 10.], [0., 1.]]);
391        let lhs = Lhs::new(&xlimits).with_rng(Xoshiro256Plus::seed_from_u64(42));
392
393        let sample1 = lhs.sample(5);
394        let sample2 = lhs.sample(5);
395        assert_abs_diff_ne!(sample1, sample2);
396    }
397
398    #[test]
399    fn test_lhs_clone_different() {
400        let xlimits = array![[-1., 1.]];
401        let rng = Xoshiro256Plus::seed_from_u64(42);
402        let lhs = Lhs::new(&xlimits)
403            .kind(LhsKind::Classic)
404            .with_rng(rng.clone());
405        let lhs1 = lhs.clone();
406        let s1 = lhs1.sample(10);
407        let lhs2 = lhs.clone();
408        let s2 = lhs2.sample(10);
409        assert_abs_diff_ne!(s1, s2);
410    }
411}
egobox_doe/lhs.rs

egobox_doe/
lhs.rs
