non_convex_opt/algorithms/multi_swarm/
particle.rs

1use crate::utils::opt_prob::{FloatNumber as FloatNum, OptProb};
2use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, OVector, U1};
3use rand::{rngs::StdRng, Rng};
4
5pub struct Particle<T, D>
6where
7    T: FloatNum + Send + Sync,
8    D: Dim + Send + Sync,
9    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<U1>,
10{
11    pub position: OVector<T, D>,
12    pub velocity: OVector<T, D>,
13    pub best_position: OVector<T, D>,
14    pub best_fitness: T,
15    pub improvement_counter: usize,
16    pub stagnation_counter: usize,
17    rng: StdRng,
18}
19
20impl<T, D> Particle<T, D>
21where
22    T: FloatNum + Send + Sync,
23    D: Dim + Send + Sync,
24    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<U1>,
25{
26    pub fn new(position: OVector<T, D>, velocity: OVector<T, D>, fitness: T, rng: StdRng) -> Self {
27        Self {
28            position: position.clone(),
29            velocity,
30            best_position: position,
31            best_fitness: fitness,
32            improvement_counter: 0,
33            stagnation_counter: 0,
34            rng,
35        }
36    }
37
38    pub fn update_velocity_and_position(
39        &mut self,
40        global_best: &OVector<T, D>,
41        w: T,
42        c1: T,
43        c2: T,
44        opt_prob: &OptProb<T, D>,
45        bounds: (T, T),
46    ) {
47        for i in 0..self.velocity.len() {
48            let r1 = T::from_f64(self.rng.random::<f64>()).unwrap();
49            let r2 = T::from_f64(self.rng.random::<f64>()).unwrap();
50
51            let cognitive = c1 * r1 * (self.best_position[i] - self.position[i]);
52            let social = c2 * r2 * (global_best[i] - self.position[i]);
53
54            // Clamp based on search space size
55            let search_space_size = bounds.1 - bounds.0;
56            let v_max = search_space_size * T::from_f64(0.2).unwrap(); // More aggressive clamping
57
58            self.velocity[i] = (w * self.velocity[i] + cognitive + social).clamp(-v_max, v_max);
59        }
60
61        // Reflective boundary
62        let new_positions: Vec<T> = self
63            .position
64            .iter()
65            .zip(self.velocity.iter())
66            .map(|(&p, &v)| {
67                let new_pos = p + v;
68                if new_pos < bounds.0 {
69                    let reflected_pos = bounds.0 + (bounds.0 - new_pos);
70                    reflected_pos.clamp(bounds.0, bounds.1)
71                } else if new_pos > bounds.1 {
72                    let reflected_pos = bounds.1 - (new_pos - bounds.1);
73                    reflected_pos.clamp(bounds.0, bounds.1)
74                } else {
75                    new_pos
76                }
77            })
78            .collect();
79
80        let final_position = OVector::<T, D>::from_vec_generic(
81            D::from_usize(new_positions.len()),
82            U1,
83            new_positions,
84        );
85        self.position = final_position;
86
87        // Only update best position if new position is better AND feasible
88        let new_fitness = opt_prob.evaluate(&self.position);
89        if new_fitness > self.best_fitness && opt_prob.is_feasible(&self.position) {
90            self.best_fitness = new_fitness;
91            self.best_position = self.position.clone();
92            self.improvement_counter += 1;
93            self.stagnation_counter = 0;
94        } else {
95            self.stagnation_counter += 1;
96        }
97    }
98
99    pub fn is_stagnated(&self, threshold: usize) -> bool {
100        self.stagnation_counter > threshold
101    }
102
103    // Iterations since last improvement
104    pub fn age(&self) -> usize {
105        self.stagnation_counter
106    }
107
108    pub fn reset_stagnation(&mut self) {
109        self.stagnation_counter = 0;
110    }
111}