rlevo_evolution/algorithms/metaheuristic/
salp.rs

1//! Salp Swarm Algorithm.
2//!
3//! The population splits into a *leader* half and a *follower* chain:
4//!
5//! - Leaders (indices `0 .. N/2`) update toward the food source `F`
6//!   (the best-so-far position), modulated by a time-decayed coefficient
7//!   `c1 = 2·exp(−(4·t / T)²)`:
8//!
9//!   `X_i ← F ± c1 · ((ub − lb)·c2 + lb)` with `c2, c3 ∈ U[0, 1]`,
10//!   where the sign follows `c3 ≥ 0.5 → +` and `c3 < 0.5 → −`.
11//!
12//! - Followers (indices `N/2 .. N`) track the chain:
13//!
14//!   `X_i ← (X_i + X_{i−1}) / 2`.
15//!
16//! The follower rule is realized as a parallel stencil — the shifted
17//! copy is `concat([last_leader, followers[..−1]])` — rather than a
18//! host-side Python-style loop. That keeps the update GPU-friendly and
19//! deterministic.
20//!
21//! # Candor
22//!
23//! Legacy comparator. The leader update is a thinly-veiled biased
24//! random walk toward the best; the follower chain is a stencil average
25//! of old positions. Neither mechanism introduces a novel search
26//! dynamic not already present in PSO or DE. Ship it for API coverage;
27//! prefer CMA-ES or LSHADE when available.
28//!
29//! # References
30//!
31//! - Mirjalili, Gandomi, Mirjalili, Saremi, Faris & Mirjalili (2017),
32//!   *Salp Swarm Algorithm*.
33
34use std::marker::PhantomData;
35
36use burn::tensor::{Distribution, Int, Tensor, TensorData, backend::Backend};
37use rand::Rng;
38use rand::RngExt;
39
40use crate::rng::{SeedPurpose, seed_stream};
41use crate::strategy::{Strategy, StrategyMetrics};
42
43/// Static configuration for [`SalpSwarm`].
44#[derive(Debug, Clone)]
45pub struct SalpConfig {
46    /// Swarm size; must be `≥ 2` so there is at least one leader and
47    /// one follower.
48    pub pop_size: usize,
49    /// Genome dimensionality.
50    pub genome_dim: usize,
51    /// Search-space bounds. The leader update pulls from a scaled
52    /// uniform over this range.
53    pub bounds: (f32, f32),
54    /// Budget pacing the `c1` decay.
55    pub max_generations: usize,
56}
57
58impl SalpConfig {
59    /// Default configuration for a given population size and genome dimensionality.
60    #[must_use]
61    pub fn default_for(pop_size: usize, genome_dim: usize) -> Self {
62        Self {
63            pop_size,
64            genome_dim,
65            bounds: (-5.12, 5.12),
66            max_generations: 500,
67        }
68    }
69}
70
71/// Generation state for [`SalpSwarm`].
72#[derive(Debug, Clone)]
73pub struct SalpState<B: Backend> {
74    /// Current positions, shape `(pop_size, D)`.
75    pub positions: Tensor<B, 2>,
76    /// Host-side fitness cache.
77    pub fitness: Vec<f32>,
78    /// Food-source (best) genome.
79    pub best_genome: Option<Tensor<B, 2>>,
80    /// Food-source fitness.
81    pub best_fitness: f32,
82    /// Generation counter.
83    pub generation: usize,
84}
85
86/// Salp Swarm Algorithm strategy.
87///
88/// # Panics
89///
90/// [`Strategy::init`] panics if `params.pop_size < 2`, since the
91/// leader/follower split requires at least one of each.
92///
93/// # Example
94///
95/// ```no_run
96/// use burn::backend::NdArray;
97/// use rlevo_evolution::algorithms::metaheuristic::salp::{SalpConfig, SalpSwarm};
98///
99/// let strategy = SalpSwarm::<NdArray>::new();
100/// let params = SalpConfig::default_for(32, 10);
101/// let _ = (strategy, params);
102/// ```
103#[derive(Debug, Clone, Copy, Default)]
104pub struct SalpSwarm<B: Backend> {
105    _backend: PhantomData<fn() -> B>,
106}
107
108impl<B: Backend> SalpSwarm<B> {
109    /// Builds a new (stateless) strategy object.
110    #[must_use]
111    pub fn new() -> Self {
112        Self {
113            _backend: PhantomData,
114        }
115    }
116}
117
118impl<B: Backend> Strategy<B> for SalpSwarm<B>
119where
120    B::Device: Clone,
121{
122    type Params = SalpConfig;
123    type State = SalpState<B>;
124    type Genome = Tensor<B, 2>;
125
126    fn init(&self, params: &SalpConfig, rng: &mut dyn Rng, device: &B::Device) -> SalpState<B> {
127        assert!(params.pop_size >= 2, "SSA requires pop_size >= 2");
128        let (lo, hi) = params.bounds;
129        B::seed(device, rng.next_u64());
130        let positions = Tensor::<B, 2>::random(
131            [params.pop_size, params.genome_dim],
132            Distribution::Uniform(f64::from(lo), f64::from(hi)),
133            device,
134        );
135        SalpState {
136            positions,
137            fitness: Vec::new(),
138            best_genome: None,
139            best_fitness: f32::INFINITY,
140            generation: 0,
141        }
142    }
143
144    fn ask(
145        &self,
146        params: &SalpConfig,
147        state: &SalpState<B>,
148        rng: &mut dyn Rng,
149        device: &B::Device,
150    ) -> (Tensor<B, 2>, SalpState<B>) {
151        if state.fitness.is_empty() {
152            return (state.positions.clone(), state.clone());
153        }
154
155        let pop_size = params.pop_size;
156        let genome_dim = params.genome_dim;
157        let n_leaders = pop_size / 2;
158        let (lo, hi) = params.bounds;
159
160        // c1 decay: 2 * exp(-(4t/T)^2)
161        #[allow(clippy::cast_precision_loss)]
162        let t = state.generation as f32;
163        #[allow(clippy::cast_precision_loss)]
164        let max_t = params.max_generations.max(1) as f32;
165        let frac = (4.0 * t / max_t).min(4.0);
166        let c1 = 2.0 * (-(frac * frac)).exp();
167
168        // Host-side sampling for c2, c3 so the leader step is
169        // backend-agnostic and reproducible under the splitmix contract.
170        let mut stream = seed_stream(rng.next_u64(), state.generation as u64, SeedPurpose::Other);
171        let mut leader_delta: Vec<f32> = Vec::with_capacity(n_leaders * genome_dim);
172        for _ in 0..n_leaders {
173            for _ in 0..genome_dim {
174                let c2: f32 = stream.random::<f32>();
175                let c3: f32 = stream.random::<f32>();
176                let scaled = (hi - lo) * c2 + lo;
177                let sign = if c3 >= 0.5 { 1.0 } else { -1.0 };
178                leader_delta.push(sign * c1 * scaled);
179            }
180        }
181
182        let best = state
183            .best_genome
184            .as_ref()
185            .expect("best_genome populated after first tell")
186            .clone()
187            .expand([n_leaders, genome_dim]);
188        let delta = Tensor::<B, 2>::from_data(
189            TensorData::new(leader_delta, [n_leaders, genome_dim]),
190            device,
191        );
192        let new_leaders = (best + delta).clamp(lo, hi);
193
194        // Follower stencil: X_i ← (X_i + X_{i-1}) / 2.
195        let followers = state
196            .positions
197            .clone()
198            .slice([n_leaders..pop_size, 0..genome_dim]);
199        // Shifted copy: follower i reads the position of the salp
200        // directly ahead, which is the last leader for the first
201        // follower and `followers[i-1]` otherwise. We build the shifted
202        // stack by gathering indices `[n_leaders-1, n_leaders, …, pop_size-2]`
203        // from the *updated* leader block + the old followers slice.
204        let joined = Tensor::cat(vec![new_leaders.clone(), followers.clone()], 0); // (pop_size, D) — leaders are already the new ones.
205        #[allow(clippy::cast_possible_wrap)]
206        let shift_idx: Vec<i64> = (0..(pop_size - n_leaders))
207            .map(|k| (n_leaders + k - 1) as i64)
208            .collect();
209        let idx = Tensor::<B, 1, Int>::from_data(
210            TensorData::new(shift_idx, [pop_size - n_leaders]),
211            device,
212        );
213        let previous = joined.clone().select(0, idx);
214        let new_followers = (followers + previous).mul_scalar(0.5).clamp(lo, hi);
215
216        let new_positions = Tensor::cat(vec![new_leaders, new_followers], 0);
217        let mut next = state.clone();
218        next.positions.clone_from(&new_positions);
219        (new_positions, next)
220    }
221
222    fn tell(
223        &self,
224        _params: &SalpConfig,
225        population: Tensor<B, 2>,
226        fitness: Tensor<B, 1>,
227        mut state: SalpState<B>,
228        _rng: &mut dyn Rng,
229    ) -> (SalpState<B>, StrategyMetrics) {
230        let fitness_host = fitness.into_data().into_vec::<f32>().unwrap_or_default();
231        state.fitness.clone_from(&fitness_host);
232        state.positions.clone_from(&population);
233        let best_idx = argmin(&fitness_host);
234        if fitness_host[best_idx] < state.best_fitness {
235            state.best_fitness = fitness_host[best_idx];
236            let device = population.device();
237            #[allow(clippy::cast_possible_wrap)]
238            let idx = Tensor::<B, 1, Int>::from_data(
239                TensorData::new(vec![best_idx as i64], [1]),
240                &device,
241            );
242            state.best_genome = Some(population.select(0, idx));
243        }
244        state.generation += 1;
245        let m =
246            StrategyMetrics::from_host_fitness(state.generation, &fitness_host, state.best_fitness);
247        state.best_fitness = m.best_fitness_ever;
248        (state, m)
249    }
250
251    fn best(&self, state: &SalpState<B>) -> Option<(Tensor<B, 2>, f32)> {
252        state
253            .best_genome
254            .as_ref()
255            .map(|g| (g.clone(), state.best_fitness))
256    }
257}
258
259fn argmin(xs: &[f32]) -> usize {
260    let mut best_idx = 0usize;
261    let mut best = f32::INFINITY;
262    for (i, &v) in xs.iter().enumerate() {
263        if v < best {
264            best = v;
265            best_idx = i;
266        }
267    }
268    best_idx
269}
270
271#[cfg(test)]
272mod tests {
273    use super::*;
274    use crate::fitness::FromFitnessEvaluable;
275    use crate::strategy::EvolutionaryHarness;
276    use burn::backend::NdArray;
277    use rlevo_core::fitness::FitnessEvaluable;
278
279    type TestBackend = NdArray;
280
281    struct Sphere;
282    struct SphereFit;
283    impl FitnessEvaluable for SphereFit {
284        type Individual = Vec<f64>;
285        type Landscape = Sphere;
286        fn evaluate(&self, x: &Self::Individual, _: &Self::Landscape) -> f64 {
287            x.iter().map(|v| v * v).sum()
288        }
289    }
290
291    #[test]
292    fn ssa_converges_on_sphere_d10() {
293        // SSA reduces strongly but, per module-level candor, is not
294        // expected to reach machine precision. A 1e-2 target in 600
295        // generations is the strong-reduction guarantee mirroring the
296        // DE premature-convergence test pattern.
297        let device = Default::default();
298        let strategy = SalpSwarm::<TestBackend>::new();
299        let params = SalpConfig::default_for(40, 10);
300        let fitness_fn = FromFitnessEvaluable::new(SphereFit, Sphere);
301        let mut harness = EvolutionaryHarness::<TestBackend, _, _>::new(
302            strategy, params, fitness_fn, 3, device, 600,
303        );
304        harness.reset();
305        while !harness.step(()).done {}
306        let best = harness.latest_metrics().unwrap().best_fitness_ever;
307        assert!(best < 1e-2, "SSA D10 best={best}");
308    }
309}
rlevo_evolution/algorithms/metaheuristic/salp.rs

rlevo_evolution/algorithms/metaheuristic/
salp.rs
