rlevo_evolution/algorithms/metaheuristic/
bat.rs

1//! Bat Algorithm.
2//!
3//! Each bat carries a position, a velocity, a frequency `f`, a loudness
4//! `A`, and a pulse rate `r`. Per generation:
5//!
6//! 1. Sample `f_i = f_min + (f_max − f_min)·β`, `β ∈ U[0, 1]`.
7//! 2. Update velocity: `v_i ← v_i + (x_i − x_best)·f_i`.
8//! 3. Propose candidate: `x'_i = x_i + v_i`. If `rand > r_i`, override
9//!    with a local walk `x'_i = x_best + ε · mean(A)`, `ε ∈ U[−1, 1]`.
10//! 4. `tell` accepts the candidate iff
11//!    `rand < A_i` **and** `f(x'_i) ≤ f(x_i)`. On acceptance:
12//!    `A_i *= α` (decay loudness), `r_i = r_{i,0}·(1 − exp(−γ·t))`
13//!    (grow pulse rate).
14//!
15//! # Candor
16//!
17//! Legacy comparator. The velocity/position update is structurally a
18//! PSO variant toward the global best; the probabilistic acceptance
19//! adds simulated-annealing-style noise. Camacho Villalón et al. (2020)
20//! discuss the lack of search mechanisms not already present in
21//! simpler algorithms. Ship it for API coverage; prefer CMA-ES or
22//! LSHADE when available.
23//!
24//! # References
25//!
26//! - Yang (2010), *A New Metaheuristic Bat-Inspired Algorithm*.
27
28use std::marker::PhantomData;
29
30use burn::tensor::{Int, Tensor, TensorData, backend::Backend};
31use rand::Rng;
32use rand::RngExt;
33
34use crate::rng::{SeedPurpose, seed_stream};
35use crate::strategy::{Strategy, StrategyMetrics};
36
37/// Static configuration for [`BatAlgorithm`].
38#[derive(Debug, Clone)]
39pub struct BatConfig {
40    /// Number of bats.
41    pub pop_size: usize,
42    /// Genome dimensionality.
43    pub genome_dim: usize,
44    /// Search-space bounds.
45    pub bounds: (f32, f32),
46    /// Minimum frequency.
47    pub f_min: f32,
48    /// Maximum frequency.
49    pub f_max: f32,
50    /// Initial loudness.
51    pub a0: f32,
52    /// Initial pulse rate.
53    pub r0: f32,
54    /// Loudness decay factor (0 < α ≤ 1). Canonical `α = 0.9`.
55    pub alpha: f32,
56    /// Pulse-rate growth factor (γ > 0). Canonical `γ = 0.9`.
57    pub gamma: f32,
58}
59
60impl BatConfig {
61    /// Default configuration for a given population size and genome dimensionality.
62    #[must_use]
63    pub fn default_for(pop_size: usize, genome_dim: usize) -> Self {
64        Self {
65            pop_size,
66            genome_dim,
67            bounds: (-5.12, 5.12),
68            f_min: 0.0,
69            f_max: 2.0,
70            a0: 1.0,
71            r0: 0.5,
72            alpha: 0.9,
73            gamma: 0.9,
74        }
75    }
76}
77
78/// Generation state for [`BatAlgorithm`].
79#[derive(Debug, Clone)]
80pub struct BatState<B: Backend> {
81    /// Current positions, shape `(pop_size, D)`.
82    pub positions: Tensor<B, 2>,
83    /// Current velocities, shape `(pop_size, D)`.
84    pub velocities: Tensor<B, 2>,
85    /// Per-bat loudness.
86    pub loudness: Vec<f32>,
87    /// Per-bat pulse rate.
88    pub pulse_rate: Vec<f32>,
89    /// Host-side fitness cache for the current positions.
90    pub fitness: Vec<f32>,
91    /// Best-so-far genome.
92    pub best_genome: Option<Tensor<B, 2>>,
93    /// Best-so-far fitness.
94    pub best_fitness: f32,
95    /// Generation counter.
96    pub generation: usize,
97    /// Per-generation "accept this candidate?" decisions recorded in
98    /// `ask` so `tell` can gate the loudness/pulse updates consistently
99    /// with the RNG draws.
100    pub pending_accept: Vec<bool>,
101}
102
103/// Bat Algorithm strategy.
104///
105/// # Example
106///
107/// ```no_run
108/// use burn::backend::Flex;
109/// use rlevo_evolution::algorithms::metaheuristic::bat::{BatAlgorithm, BatConfig};
110///
111/// let strategy = BatAlgorithm::<Flex>::new();
112/// let params = BatConfig::default_for(32, 10);
113/// let _ = (strategy, params);
114/// ```
115#[derive(Debug, Clone, Copy, Default)]
116pub struct BatAlgorithm<B: Backend> {
117    _backend: PhantomData<fn() -> B>,
118}
119
120impl<B: Backend> BatAlgorithm<B> {
121    /// Builds a new (stateless) strategy object.
122    #[must_use]
123    pub fn new() -> Self {
124        Self {
125            _backend: PhantomData,
126        }
127    }
128}
129
130impl<B: Backend> Strategy<B> for BatAlgorithm<B>
131where
132    B::Device: Clone,
133{
134    type Params = BatConfig;
135    type State = BatState<B>;
136    type Genome = Tensor<B, 2>;
137
138    /// Build the initial colony by host-sampling `pop_size` positions
139    /// uniformly in `[bounds.lo, bounds.hi]`.
140    ///
141    /// Velocities are zeroed, loudness is set to `params.a0`, pulse rate
142    /// to `params.r0`, and `fitness` left empty so that the first
143    /// [`ask`](Strategy::ask) → [`tell`](Strategy::tell) pair initialises
144    /// those fields before any acceptance logic runs.  Positions are drawn
145    /// from a deterministic [`seed_stream`]; the process-wide Flex RNG is
146    /// never touched.
147    fn init(&self, params: &BatConfig, rng: &mut dyn Rng, device: &<B as burn::tensor::backend::BackendTypes>::Device) -> BatState<B> {
148        let (lo, hi) = params.bounds;
149        // Sample initial positions on the host from a deterministic
150        // `seed_stream`, mirroring `ask`/`tell`. The Flex backend's
151        // `Tensor::random` draws from a process-wide RNG mutex; under the
152        // parallel test runner those draws interleave with sibling tests,
153        // so `B::seed` + `Tensor::random` is NOT reproducible across
154        // thread schedules. Host sampling keeps initialisation bit-stable
155        // regardless of core count or test ordering.
156        let pop = params.pop_size;
157        let genome_dim = params.genome_dim;
158        let mut stream = seed_stream(rng.next_u64(), 0, SeedPurpose::Init);
159        let mut position_rows = Vec::with_capacity(pop * genome_dim);
160        for _ in 0..pop * genome_dim {
161            position_rows.push(lo + (hi - lo) * stream.random::<f32>());
162        }
163        let positions =
164            Tensor::<B, 2>::from_data(TensorData::new(position_rows, [pop, genome_dim]), device);
165        let velocities = Tensor::<B, 2>::zeros([params.pop_size, params.genome_dim], device);
166        BatState {
167            positions,
168            velocities,
169            loudness: vec![params.a0; params.pop_size],
170            pulse_rate: vec![params.r0; params.pop_size],
171            fitness: Vec::new(),
172            best_genome: None,
173            best_fitness: f32::INFINITY,
174            generation: 0,
175            pending_accept: Vec::new(),
176        }
177    }
178
179    /// Propose candidate positions for the current generation.
180    ///
181    /// On the first call (`state.fitness` is empty) returns the initial
182    /// positions unchanged so the caller can evaluate generation zero.
183    ///
184    /// On subsequent calls the update proceeds in three host/device steps:
185    ///
186    /// 1. **Frequency** — sample `f_i = f_min + (f_max − f_min)·β_i`,
187    ///    `β_i ∈ U[0,1]`.
188    /// 2. **Global move** — `v_i ← v_i + (x_i − x_best)·f_i`,
189    ///    `x'_i = x_i + v_i`.
190    /// 3. **Local walk** (when `rand > r_i`) — override with
191    ///    `x'_i = x_best + ε·mean(A)`, `ε ∈ U[−1,1]`.
192    ///
193    /// All random draws are host-sampled through [`seed_stream`] for
194    /// bit-stable reproduction across thread schedules.  The
195    /// per-bat acceptance decisions (`pending_accept`) are recorded in the
196    /// returned state and consumed by [`tell`](Strategy::tell).
197    ///
198    /// # Panics
199    ///
200    /// Panics if called when `state.best_genome` is `None` after the first
201    /// generation has been evaluated (i.e. if `state.fitness` is non-empty
202    /// but `state.best_genome` was not set by a preceding `tell` call).
203    fn ask(
204        &self,
205        params: &BatConfig,
206        state: &BatState<B>,
207        rng: &mut dyn Rng,
208        device: &<B as burn::tensor::backend::BackendTypes>::Device,
209    ) -> (Tensor<B, 2>, BatState<B>) {
210        if state.fitness.is_empty() {
211            // Evaluate the initial colony first; the velocity update is
212            // only defined once a best exists.
213            return (state.positions.clone(), state.clone());
214        }
215
216        let pop = params.pop_size;
217        let genome_dim = params.genome_dim;
218        let (lo, hi) = params.bounds;
219
220        // Host-side sampling for β, pulse check, acceptance draw, and
221        // local-walk ε. Keeping these on host preserves bit-parity
222        // across backends (Mantegna / wgpu normal RNG has documented
223        // fp drift under FMA reordering; BA draws are mostly uniform).
224        let mut stream = seed_stream(rng.next_u64(), state.generation as u64, SeedPurpose::Other);
225
226        let mut betas = Vec::with_capacity(pop);
227        let mut use_local = Vec::with_capacity(pop);
228        let mut accept_draw = Vec::with_capacity(pop);
229        let mut epsilon_rows = Vec::with_capacity(pop * genome_dim);
230        for i in 0..pop {
231            betas.push(stream.random::<f32>());
232            use_local.push(stream.random::<f32>() > state.pulse_rate[i]);
233            accept_draw.push(stream.random::<f32>());
234            for _ in 0..genome_dim {
235                epsilon_rows.push(2.0 * stream.random::<f32>() - 1.0);
236            }
237        }
238
239        // Mean loudness across the colony — used by the local-walk
240        // step to scale its ε perturbation.
241        let mean_loudness: f32 = {
242            let s: f32 = state.loudness.iter().sum();
243            #[allow(clippy::cast_precision_loss)]
244            {
245                s / pop as f32
246            }
247        };
248
249        let best = state
250            .best_genome
251            .as_ref()
252            .expect("best populated after first tell")
253            .clone()
254            .expand([pop, genome_dim]);
255
256        // Frequency: f_i = f_min + (f_max - f_min) · β_i  → shape (pop, 1) → (pop, D).
257        let f_vec: Vec<f32> = betas
258            .iter()
259            .map(|b| params.f_min + (params.f_max - params.f_min) * b)
260            .collect();
261        let f_mat = Tensor::<B, 1>::from_data(TensorData::new(f_vec, [pop]), device)
262            .unsqueeze_dim::<2>(1)
263            .expand([pop, genome_dim]);
264
265        let new_velocities =
266            state.velocities.clone() + (state.positions.clone() - best.clone()).mul(f_mat);
267        let global_move = state.positions.clone() + new_velocities.clone();
268        // Local walk: x_best + ε · mean(A).
269        let eps =
270            Tensor::<B, 2>::from_data(TensorData::new(epsilon_rows, [pop, genome_dim]), device);
271        let local_move = best + eps.mul_scalar(mean_loudness);
272
273        #[allow(clippy::cast_possible_wrap)]
274        let mask = Tensor::<B, 1, Int>::from_data(
275            TensorData::new(
276                use_local.iter().map(|&b| i64::from(b)).collect::<Vec<_>>(),
277                [pop],
278            ),
279            device,
280        )
281        .equal_elem(1)
282        .unsqueeze_dim::<2>(1)
283        .expand([pop, genome_dim]);
284        let candidates = global_move.mask_where(mask, local_move).clamp(lo, hi);
285
286        // Defer acceptance decisions to tell; record the random draws.
287        let mut next = state.clone();
288        next.velocities = new_velocities;
289        next.pending_accept = accept_draw
290            .iter()
291            .zip(state.loudness.iter())
292            .map(|(&draw, &a)| draw < a)
293            .collect();
294        (candidates, next)
295    }
296
297    /// Ingest candidate fitness values, apply the acceptance gate, and
298    /// advance the generation counter.
299    ///
300    /// On the first call (generation zero bootstrap) all candidates are
301    /// unconditionally accepted and loudness/pulse-rate updates are
302    /// skipped.
303    ///
304    /// On subsequent calls candidate `i` replaces position `i` iff
305    /// `pending_accept[i]` (drawn in [`ask`](Strategy::ask)) **and**
306    /// `fitness[i] ≤ state.fitness[i]`.  On acceptance, loudness decays
307    /// (`A_i *= α`) and pulse rate grows
308    /// (`r_i = r₀·(1 − exp(−γ·t))`).
309    fn tell(
310        &self,
311        params: &BatConfig,
312        candidates: Tensor<B, 2>,
313        fitness: Tensor<B, 1>,
314        mut state: BatState<B>,
315        _rng: &mut dyn Rng,
316    ) -> (BatState<B>, StrategyMetrics) {
317        let fitness_host = fitness.into_data().into_vec::<f32>().unwrap_or_default();
318        let device = candidates.device();
319        let pop = params.pop_size;
320        let genome_dim = params.genome_dim;
321
322        if state.fitness.is_empty() {
323            state.fitness.clone_from(&fitness_host);
324            let best_idx = argmin(&fitness_host);
325            state.best_fitness = fitness_host[best_idx];
326            #[allow(clippy::cast_possible_wrap)]
327            let idx = Tensor::<B, 1, Int>::from_data(
328                TensorData::new(vec![best_idx as i64], [1]),
329                &device,
330            );
331            state.best_genome = Some(candidates.clone().select(0, idx));
332            state.positions = candidates;
333            state.generation += 1;
334            let m = StrategyMetrics::from_host_fitness(
335                state.generation,
336                &fitness_host,
337                state.best_fitness,
338            );
339            state.best_fitness = m.best_fitness_ever;
340            return (state, m);
341        }
342
343        // Acceptance: accept candidate `i` iff `pending_accept[i]` AND
344        // candidate's fitness is no worse than current.
345        #[allow(clippy::cast_possible_wrap)]
346        let mut rs: Vec<i64> = (0..pop).map(|i| i as i64).collect();
347        let mut new_fitness = state.fitness.clone();
348        #[allow(clippy::cast_precision_loss)]
349        let t = state.generation as f32;
350        for i in 0..pop {
351            let accept_gate = state.pending_accept.get(i).copied().unwrap_or(false);
352            let improves = fitness_host[i] <= state.fitness[i];
353            if accept_gate && improves {
354                #[allow(clippy::cast_possible_wrap)]
355                {
356                    rs[i] = (pop + i) as i64;
357                }
358                new_fitness[i] = fitness_host[i];
359                state.loudness[i] *= params.alpha;
360                state.pulse_rate[i] = params.r0 * (1.0 - (-params.gamma * t).exp());
361            }
362        }
363        let stacked = Tensor::cat(vec![state.positions.clone(), candidates], 0);
364        let idx = Tensor::<B, 1, Int>::from_data(TensorData::new(rs, [pop]), &device);
365        state.positions = stacked.select(0, idx);
366        state.fitness = new_fitness;
367
368        // Refresh global best.
369        let best_idx = argmin(&state.fitness);
370        if state.fitness[best_idx] < state.best_fitness {
371            state.best_fitness = state.fitness[best_idx];
372            #[allow(clippy::cast_possible_wrap)]
373            let idx = Tensor::<B, 1, Int>::from_data(
374                TensorData::new(vec![best_idx as i64], [1]),
375                &device,
376            );
377            state.best_genome = Some(state.positions.clone().select(0, idx));
378        }
379
380        state.generation += 1;
381        let m =
382            StrategyMetrics::from_host_fitness(state.generation, &fitness_host, state.best_fitness);
383        state.best_fitness = m.best_fitness_ever;
384        let _ = genome_dim;
385        (state, m)
386    }
387
388    /// Returns the best-so-far `(genome, fitness)` pair, or `None` before
389    /// the first [`tell`](Strategy::tell) call.
390    fn best(&self, state: &BatState<B>) -> Option<(Tensor<B, 2>, f32)> {
391        state
392            .best_genome
393            .as_ref()
394            .map(|g| (g.clone(), state.best_fitness))
395    }
396}
397
398fn argmin(xs: &[f32]) -> usize {
399    let mut best_idx = 0usize;
400    let mut best = f32::INFINITY;
401    for (i, &v) in xs.iter().enumerate() {
402        if v < best {
403            best = v;
404            best_idx = i;
405        }
406    }
407    best_idx
408}
409
410#[cfg(test)]
411mod tests {
412    use super::*;
413    use crate::fitness::FromFitnessEvaluable;
414    use crate::strategy::EvolutionaryHarness;
415    use burn::backend::Flex;
416    use rlevo_core::fitness::FitnessEvaluable;
417
418    type TestBackend = Flex;
419
420    struct Sphere;
421    struct SphereFit;
422    impl FitnessEvaluable for SphereFit {
423        type Individual = Vec<f64>;
424        type Landscape = Sphere;
425        fn evaluate(&self, x: &Self::Individual, _: &Self::Landscape) -> f64 {
426            x.iter().map(|v| v * v).sum()
427        }
428    }
429
430    #[test]
431    fn bat_converges_on_sphere_d10() {
432        // Bat is a "legacy comparator" per the module-level candor
433        // note. We require strong reduction from the random baseline,
434        // not machine precision — the probabilistic acceptance gate
435        // (A_i decay) throttles late-stage progress. Threshold 0.1 on
436        // Sphere-D10 in 800 generations is still orders of magnitude
437        // below the uniform-random baseline (≈ 87).
438        let device = Default::default();
439        let strategy = BatAlgorithm::<TestBackend>::new();
440        let params = BatConfig::default_for(40, 10);
441        let fitness_fn = FromFitnessEvaluable::new(SphereFit, Sphere);
442        let mut harness = EvolutionaryHarness::<TestBackend, _, _>::new(
443            strategy, params, fitness_fn, 23, device, 800,
444        );
445        harness.reset();
446        while !harness.step(()).done {}
447        let best = harness.latest_metrics().unwrap().best_fitness_ever;
448        assert!(best < 0.1, "Bat D10 best={best}");
449    }
450}
rlevo_evolution/algorithms/metaheuristic/bat.rs

rlevo_evolution/algorithms/metaheuristic/
bat.rs
