samyama_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)
39            .map(|_| Array1::zeros(dim))
40            .collect();
41
42        let mut history = Vec::with_capacity(self.config.max_iterations);
43        
44        let mut best_vars = population[0].variables.clone();
45        let mut best_fitness = population[0].fitness;
46
47        for iter in 0..self.config.max_iterations {
48            // Update fitness and find best/worst
49            let mut worst_fitness = population[0].fitness;
50            for ind in &population {
51                if ind.fitness < best_fitness {
52                    best_fitness = ind.fitness;
53                    best_vars = ind.variables.clone();
54                }
55                if ind.fitness > worst_fitness {
56                    worst_fitness = ind.fitness;
57                }
58            }
59            history.push(best_fitness);
60
61            // 2. Calculate Masses
62            let mut masses = vec![0.0; n];
63            let mut total_m = 0.0;
64            
65            // Normalize fitness to calculate mass
66            if (worst_fitness - best_fitness).abs() > 1e-9 {
67                for i in 0..n {
68                    masses[i] = (population[i].fitness - worst_fitness) / (best_fitness - worst_fitness);
69                    total_m += masses[i];
70                }
71                for i in 0..n {
72                    masses[i] /= total_m;
73                }
74            } else {
75                for i in 0..n {
76                    masses[i] = 1.0 / (n as f64);
77                }
78            }
79
80            // 3. Calculate Force and Acceleration
81            let g = self.g0 * (-(self.alpha * (iter as f64) / (self.config.max_iterations as f64))).exp();
82            let mut accelerations: Vec<Array1<f64>> = (0..n).map(|_| Array1::zeros(dim)).collect();
83
84            for i in 0..n {
85                for j in 0..n {
86                    if i == j { continue; }
87                    
88                    let mut dist_sq = 0.0;
89                    for k in 0..dim {
90                        let diff = population[j].variables[k] - population[i].variables[k];
91                        dist_sq += diff * diff;
92                    }
93                    let dist = dist_sq.sqrt() + 1e-6; // avoid zero
94
95                    for k in 0..dim {
96                        let force = g * (masses[i] * masses[j] / dist) * (population[j].variables[k] - population[i].variables[k]);
97                        accelerations[i][k] += rng.gen::<f64>() * force;
98                    }
99                }
100            }
101
102            // 4. Update Velocity and Position
103            for i in 0..n {
104                for k in 0..dim {
105                    velocities[i][k] = rng.gen::<f64>() * velocities[i][k] + accelerations[i][k];
106                    population[i].variables[k] = (population[i].variables[k] + velocities[i][k]).clamp(lower[k], upper[k]);
107                }
108                population[i].fitness = problem.fitness(&population[i].variables);
109            }
110        }
111
112        OptimizationResult {
113            best_variables: best_vars,
114            best_fitness,
115            history,
116        }
117    }
118}
samyama_optimization/algorithms/gsa.rs

samyama_optimization/algorithms/
gsa.rs
