graphmind_optimization/algorithms/
gsa.rs

1use crate::common::{Individual, OptimizationResult, Problem, SolverConfig};
2use ndarray::Array1;
3use rand::prelude::*;
4
5pub struct GSASolver {
6    pub config: SolverConfig,
7    pub g0: f64,    // Initial gravitational constant (100.0)
8    pub alpha: f64, // Damping ratio (20.0)
9}
10
11impl GSASolver {
12    pub fn new(config: SolverConfig) -> Self {
13        Self {
14            config,
15            g0: 100.0,
16            alpha: 20.0,
17        }
18    }
19
20    pub fn solve<P: Problem>(&self, problem: &P) -> OptimizationResult {
21        let mut rng = thread_rng();
22        let dim = problem.dim();
23        let (lower, upper) = problem.bounds();
24        let n = self.config.population_size;
25
26        // 1. Initialize Agents
27        let mut population: Vec<Individual> = (0..n)
28            .map(|_| {
29                let mut vars = Array1::zeros(dim);
30                for i in 0..dim {
31                    vars[i] = rng.gen_range(lower[i]..upper[i]);
32                }
33                let fitness = problem.fitness(&vars);
34                Individual::new(vars, fitness)
35            })
36            .collect();
37
38        let mut velocities: Vec<Array1<f64>> = (0..n).map(|_| Array1::zeros(dim)).collect();
39
40        let mut history = Vec::with_capacity(self.config.max_iterations);
41
42        let mut best_vars = population[0].variables.clone();
43        let mut best_fitness = population[0].fitness;
44
45        for iter in 0..self.config.max_iterations {
46            // Update fitness and find best/worst
47            let mut worst_fitness = population[0].fitness;
48            for ind in &population {
49                if ind.fitness < best_fitness {
50                    best_fitness = ind.fitness;
51                    best_vars = ind.variables.clone();
52                }
53                if ind.fitness > worst_fitness {
54                    worst_fitness = ind.fitness;
55                }
56            }
57            history.push(best_fitness);
58
59            // 2. Calculate Masses
60            let mut masses = vec![0.0; n];
61            let mut total_m = 0.0;
62
63            // Normalize fitness to calculate mass
64            if (worst_fitness - best_fitness).abs() > 1e-9 {
65                for (mass, ind) in masses.iter_mut().zip(population.iter()) {
66                    *mass = (ind.fitness - worst_fitness) / (best_fitness - worst_fitness);
67                    total_m += *mass;
68                }
69                for mass in masses.iter_mut() {
70                    *mass /= total_m;
71                }
72            } else {
73                for mass in masses.iter_mut() {
74                    *mass = 1.0 / (n as f64);
75                }
76            }
77
78            // 3. Calculate Force and Acceleration
79            let g = self.g0
80                * (-(self.alpha * (iter as f64) / (self.config.max_iterations as f64))).exp();
81            let mut accelerations: Vec<Array1<f64>> = (0..n).map(|_| Array1::zeros(dim)).collect();
82
83            for i in 0..n {
84                for j in 0..n {
85                    if i == j {
86                        continue;
87                    }
88
89                    let mut dist_sq = 0.0;
90                    for k in 0..dim {
91                        let diff = population[j].variables[k] - population[i].variables[k];
92                        dist_sq += diff * diff;
93                    }
94                    let dist = dist_sq.sqrt() + 1e-6; // avoid zero
95
96                    for (k, acc) in accelerations[i].iter_mut().enumerate().take(dim) {
97                        let force = g
98                            * (masses[i] * masses[j] / dist)
99                            * (population[j].variables[k] - population[i].variables[k]);
100                        *acc += rng.gen::<f64>() * force;
101                    }
102                }
103            }
104
105            // 4. Update Velocity and Position
106            for i in 0..n {
107                for k in 0..dim {
108                    velocities[i][k] = rng.gen::<f64>() * velocities[i][k] + accelerations[i][k];
109                    population[i].variables[k] =
110                        (population[i].variables[k] + velocities[i][k]).clamp(lower[k], upper[k]);
111                }
112                population[i].fitness = problem.fitness(&population[i].variables);
113            }
114        }
115
116        OptimizationResult {
117            best_variables: best_vars,
118            best_fitness,
119            history,
120        }
121    }
122}
graphmind_optimization/algorithms/gsa.rs

graphmind_optimization/algorithms/
gsa.rs
