pmcore/routines/initialization/
latin.rs

1use anyhow::Result;
2use faer::Mat;
3use rand::prelude::*;
4use rand::rngs::StdRng;
5use rand::Rng;
6
7use crate::prelude::Parameters;
8use crate::structs::theta::Theta;
9
10/// Generates an instance of [Theta] using Latin Hypercube Sampling.
11///
12/// # Arguments
13///
14/// * `parameters` - The [Parameters] struct, which contains the parameters to be sampled.
15/// * `points` - The number of points to generate, i.e. the number of rows in the matrix.
16/// * `seed` - The seed for the Sobol sequence generator.
17///
18/// # Returns
19///
20/// [Theta], a structure that holds the support point matrix
21///
22pub fn generate(parameters: &Parameters, points: usize, seed: usize) -> Result<Theta> {
23    let params: Vec<(String, f64, f64, bool)> = parameters
24        .iter()
25        .map(|p| (p.name.clone(), p.lower, p.upper, p.fixed))
26        .collect();
27
28    // Random parameters are sampled from the Sobol sequence
29    let random_params: Vec<(String, f64, f64)> = params
30        .iter()
31        .filter(|(_, _, _, fixed)| !fixed)
32        .map(|(name, lower, upper, _)| (name.clone(), *lower, *upper))
33        .collect();
34
35    // Initialize random number generator with the provided seed
36    let mut rng = StdRng::seed_from_u64(seed as u64);
37
38    // Create and shuffle intervals for each parameter
39    let mut intervals = Vec::new();
40    for _ in 0..random_params.len() {
41        let mut param_intervals: Vec<f64> = (0..points).map(|i| i as f64).collect();
42        param_intervals.shuffle(&mut rng);
43        intervals.push(param_intervals);
44    }
45
46    let rand_matrix = Mat::from_fn(points, random_params.len(), |i, j| {
47        // Get the interval for this parameter and point
48        let interval = intervals[j][i];
49        let random_offset = rng.random::<f64>();
50        // Calculate normalized value in [0,1]
51        let unscaled = (interval + random_offset) / points as f64;
52        // Scale to parameter range
53        let (_name, lower, upper) = random_params.get(j).unwrap(); // Fixed: use j instead of i
54        lower + unscaled * (upper - lower)
55    });
56
57    // Fixed parameters are initialized to the middle of their range
58    let fixed_params: Vec<(String, f64)> = params
59        .iter()
60        .filter(|(_, _, _, fixed)| *fixed)
61        .map(|(name, lower, upper, _)| (name.clone(), (upper - lower) / 2.0))
62        .collect();
63
64    let theta = Theta::from_parts(rand_matrix, random_params, fixed_params);
65
66    Ok(theta)
67}
pmcore/routines/initialization/latin.rs

pmcore/routines/initialization/
latin.rs
