type HillClimbOptimizer = fn(fn(&[f64]) -> f64, &[f64], &[(f64, f64)], f64, usize, u64) -> Vec<f64>;
pub fn simulated_annealing(
objective: fn(&[f64]) -> f64,
x0: &[f64],
bounds: &[(f64, f64)],
t_initial: f64,
t_min: f64,
cooling_rate: f64,
steps_per_temp: usize,
seed: u64,
) -> Vec<f64> {
let mut rng = SimpleRng::new(seed);
let mut x = x0.to_vec();
let mut best = x.clone();
let mut e = objective(&x);
let mut e_best = e;
let mut temp = t_initial;
while temp > t_min {
for _ in 0..steps_per_temp {
let mut candidate = x.clone();
let idx = rng.next_usize(x.len());
let range = bounds[idx].1 - bounds[idx].0;
candidate[idx] += (rng.next_f64() - 0.5) * range * 0.1;
candidate[idx] = candidate[idx].clamp(bounds[idx].0, bounds[idx].1);
let e_new = objective(&candidate);
let delta = e_new - e;
if delta < 0.0 || rng.next_f64() < (-delta / temp).exp() {
x = candidate;
e = e_new;
if e < e_best {
best = x.clone();
e_best = e;
}
}
}
temp *= cooling_rate;
}
best
}
pub fn particle_swarm(
objective: fn(&[f64]) -> f64,
bounds: &[(f64, f64)],
n_particles: usize,
max_iter: usize,
w: f64,
c1: f64,
c2: f64,
seed: u64,
) -> Vec<f64> {
let n = bounds.len();
let mut rng = SimpleRng::new(seed);
let mut positions: Vec<Vec<f64>> = (0..n_particles)
.map(|_| {
(0..n)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect()
})
.collect();
let mut velocities: Vec<Vec<f64>> = (0..n_particles).map(|_| vec![0.0; n]).collect();
let mut p_best = positions.clone();
let mut p_best_val: Vec<f64> = positions.iter().map(|p| objective(p)).collect();
let g_best_idx = p_best_val
.iter()
.enumerate()
.min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
.unwrap()
.0;
let mut g_best = positions[g_best_idx].clone();
let mut g_best_val = p_best_val[g_best_idx];
for _ in 0..max_iter {
for i in 0..n_particles {
for j in 0..n {
let r1 = rng.next_f64();
let r2 = rng.next_f64();
velocities[i][j] = w * velocities[i][j]
+ c1 * r1 * (p_best[i][j] - positions[i][j])
+ c2 * r2 * (g_best[j] - positions[i][j]);
positions[i][j] =
(positions[i][j] + velocities[i][j]).clamp(bounds[j].0, bounds[j].1);
}
let val = objective(&positions[i]);
if val < p_best_val[i] {
p_best_val[i] = val;
p_best[i] = positions[i].clone();
if val < g_best_val {
g_best_val = val;
g_best = positions[i].clone();
}
}
}
}
g_best
}
pub fn nelder_mead(
f: fn(&[f64]) -> f64,
x0: &[f64],
step: f64,
max_iter: usize,
tol: f64,
) -> Vec<f64> {
let n = x0.len();
let mut simplex: Vec<Vec<f64>> = Vec::with_capacity(n + 1);
simplex.push(x0.to_vec());
for i in 0..n {
let mut v = x0.to_vec();
v[i] += step;
simplex.push(v);
}
for _ in 0..max_iter {
let mut order: Vec<usize> = (0..=n).collect();
order.sort_by(|&a, &b| {
f(&simplex[a])
.partial_cmp(&f(&simplex[b]))
.unwrap_or(std::cmp::Ordering::Equal)
});
let best_val = f(&simplex[order[0]]);
let worst_val = f(&simplex[order[n]]);
if worst_val - best_val < tol {
break;
}
let centroid: Vec<f64> = (0..n)
.map(|j| order[..n].iter().map(|&i| simplex[i][j]).sum::<f64>() / n as f64)
.collect();
let reflected: Vec<f64> = (0..n)
.map(|j| 2.0 * centroid[j] - simplex[order[n]][j])
.collect();
let fr = f(&reflected);
if fr < f(&simplex[order[n - 1]]) && fr >= best_val {
simplex[order[n]] = reflected;
} else if fr < best_val {
let expanded: Vec<f64> = (0..n)
.map(|j| 3.0 * centroid[j] - 2.0 * simplex[order[n]][j])
.collect();
if f(&expanded) < fr {
simplex[order[n]] = expanded;
} else {
simplex[order[n]] = reflected;
}
} else {
let contracted: Vec<f64> = (0..n)
.map(|j| 0.5 * (centroid[j] + simplex[order[n]][j]))
.collect();
if f(&contracted) < worst_val {
simplex[order[n]] = contracted;
} else {
for i in 1..=n {
let best = simplex[order[0]].clone();
for (sij, &bj) in simplex[order[i]].iter_mut().zip(&best) {
*sij = 0.5 * (*sij + bj);
}
}
}
}
}
let mut order: Vec<usize> = (0..=n).collect();
order.sort_by(|&a, &b| {
f(&simplex[a])
.partial_cmp(&f(&simplex[b]))
.unwrap_or(std::cmp::Ordering::Equal)
});
simplex[order[0]].clone()
}
struct SimpleRng {
state: u64,
}
impl SimpleRng {
fn new(seed: u64) -> Self {
Self { state: seed }
}
fn next_u64(&mut self) -> u64 {
self.state = self
.state
.wrapping_mul(6364136223846793005)
.wrapping_add(1442695040888963407);
self.state
}
fn next_f64(&mut self) -> f64 {
(self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
}
fn next_usize(&mut self, max: usize) -> usize {
(self.next_f64() * max as f64) as usize % max
}
}
pub fn tabu_search(
objective: fn(&[f64]) -> f64,
x0: &[f64],
bounds: &[(f64, f64)],
tabu_size: usize,
n_neighbors: usize,
max_iter: usize,
seed: u64,
) -> Vec<f64> {
let n = x0.len();
let mut rng = SimpleRng::new(seed);
let mut best = x0.to_vec();
let mut best_val = objective(&best);
let mut current = best.clone();
let mut tabu_list: Vec<Vec<f64>> = Vec::new();
for _ in 0..max_iter {
let mut best_neighbor = current.clone();
let mut best_neighbor_val = f64::INFINITY;
for _ in 0..n_neighbors {
let mut neighbor = current.clone();
let idx = rng.next_usize(n);
let range = bounds[idx].1 - bounds[idx].0;
neighbor[idx] = (neighbor[idx] + (rng.next_f64() - 0.5) * range * 0.2)
.clamp(bounds[idx].0, bounds[idx].1);
let is_tabu = tabu_list.iter().any(|t| {
t.iter()
.zip(neighbor.iter())
.all(|(a, b)| (a - b).abs() < 1e-10)
});
if is_tabu {
continue;
}
let val = objective(&neighbor);
if val < best_neighbor_val {
best_neighbor = neighbor;
best_neighbor_val = val;
}
}
current = best_neighbor.clone();
tabu_list.push(best_neighbor.clone());
if tabu_list.len() > tabu_size {
tabu_list.remove(0);
}
if best_neighbor_val < best_val {
best = best_neighbor;
best_val = best_neighbor_val;
}
}
best
}
pub fn firefly_algorithm(
objective: fn(&[f64]) -> f64,
bounds: &[(f64, f64)],
pop_size: usize,
max_iter: usize,
alpha: f64,
beta0: f64,
gamma: f64,
seed: u64,
) -> Vec<f64> {
let n = bounds.len();
let mut rng = SimpleRng::new(seed);
let mut pop: Vec<Vec<f64>> = (0..pop_size)
.map(|_| {
(0..n)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect()
})
.collect();
let mut fitness: Vec<f64> = pop.iter().map(|p| objective(p)).collect();
for _ in 0..max_iter {
for i in 0..pop_size {
for j in 0..pop_size {
if fitness[j] < fitness[i] {
let dist_sq: f64 = pop[i]
.iter()
.zip(pop[j].iter())
.map(|(a, b)| (a - b) * (a - b))
.sum();
let beta = beta0 * (-gamma * dist_sq).exp();
for (k, bk) in bounds.iter().enumerate() {
pop[i][k] +=
beta * (pop[j][k] - pop[i][k]) + alpha * (rng.next_f64() - 0.5);
pop[i][k] = pop[i][k].clamp(bk.0, bk.1);
}
fitness[i] = objective(&pop[i]);
}
}
}
}
let best = fitness
.iter()
.enumerate()
.min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
.unwrap()
.0;
pop[best].clone()
}
pub fn grey_wolf_optimizer(
objective: fn(&[f64]) -> f64,
bounds: &[(f64, f64)],
pop_size: usize,
max_iter: usize,
seed: u64,
) -> Vec<f64> {
let n = bounds.len();
let mut rng = SimpleRng::new(seed);
let mut pop: Vec<Vec<f64>> = (0..pop_size)
.map(|_| {
(0..n)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect()
})
.collect();
let mut fitness: Vec<f64> = pop.iter().map(|p| objective(p)).collect();
for iter in 0..max_iter {
let mut sorted: Vec<usize> = (0..pop_size).collect();
sorted.sort_by(|&a, &b| {
fitness[a]
.partial_cmp(&fitness[b])
.unwrap_or(std::cmp::Ordering::Equal)
});
let alpha_w = pop[sorted[0]].clone();
let beta_w = pop[sorted[1.min(pop_size - 1)]].clone();
let delta_w = pop[sorted[2.min(pop_size - 1)]].clone();
let a_coeff = 2.0 * (1.0 - iter as f64 / max_iter as f64);
for i in 0..pop_size {
for j in 0..n {
let r1 = rng.next_f64();
let r2 = rng.next_f64();
let a1 = 2.0 * a_coeff * r1 - a_coeff;
let c1 = 2.0 * r2;
let d_alpha = (c1 * alpha_w[j] - pop[i][j]).abs();
let x1 = alpha_w[j] - a1 * d_alpha;
let r1 = rng.next_f64();
let r2 = rng.next_f64();
let a2 = 2.0 * a_coeff * r1 - a_coeff;
let c2 = 2.0 * r2;
let d_beta = (c2 * beta_w[j] - pop[i][j]).abs();
let x2 = beta_w[j] - a2 * d_beta;
let r1 = rng.next_f64();
let r2 = rng.next_f64();
let a3 = 2.0 * a_coeff * r1 - a_coeff;
let c3 = 2.0 * r2;
let d_delta = (c3 * delta_w[j] - pop[i][j]).abs();
let x3 = delta_w[j] - a3 * d_delta;
pop[i][j] = ((x1 + x2 + x3) / 3.0).clamp(bounds[j].0, bounds[j].1);
}
fitness[i] = objective(&pop[i]);
}
}
let best = fitness
.iter()
.enumerate()
.min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
.unwrap()
.0;
pop[best].clone()
}
pub fn hill_climbing(
objective: fn(&[f64]) -> f64,
x0: &[f64],
bounds: &[(f64, f64)],
step_size: f64,
max_iter: usize,
seed: u64,
) -> Vec<f64> {
let n = x0.len();
let mut rng = SimpleRng::new(seed);
let mut x = x0.to_vec();
let mut fx = objective(&x);
for _ in 0..max_iter {
let mut candidate = x.clone();
let idx = rng.next_usize(n);
let range = bounds[idx].1 - bounds[idx].0;
candidate[idx] = (candidate[idx] + (rng.next_f64() - 0.5) * range * step_size)
.clamp(bounds[idx].0, bounds[idx].1);
let fc = objective(&candidate);
if fc < fx {
x = candidate;
fx = fc;
}
}
x
}
pub fn random_search(
objective: fn(&[f64]) -> f64,
bounds: &[(f64, f64)],
n_samples: usize,
seed: u64,
) -> Vec<f64> {
let n = bounds.len();
let mut rng = SimpleRng::new(seed);
let mut best: Vec<f64> = (0..n)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect();
let mut best_val = objective(&best);
for _ in 1..n_samples {
let candidate: Vec<f64> = (0..n)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect();
let val = objective(&candidate);
if val < best_val {
best = candidate;
best_val = val;
}
}
best
}
pub fn multistart(
optimizer: HillClimbOptimizer,
objective: fn(&[f64]) -> f64,
bounds: &[(f64, f64)],
step_size: f64,
n_starts: usize,
iter_per_start: usize,
seed: u64,
) -> Vec<f64> {
let n = bounds.len();
let mut rng = SimpleRng::new(seed);
let mut global_best: Vec<f64> = Vec::new();
let mut global_best_val = f64::INFINITY;
for _ in 0..n_starts {
let x0: Vec<f64> = (0..n)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect();
let s = rng.next_u64();
let result = optimizer(objective, &x0, bounds, step_size, iter_per_start, s);
let val = objective(&result);
if val < global_best_val {
global_best = result;
global_best_val = val;
}
}
global_best
}
pub fn whale_optimization(
objective: fn(&[f64]) -> f64,
bounds: &[(f64, f64)],
pop_size: usize,
max_iter: usize,
seed: u64,
) -> Vec<f64> {
let dim = bounds.len();
let mut rng = SimpleRng::new(seed);
let mut pop: Vec<Vec<f64>> = (0..pop_size)
.map(|_| {
(0..dim)
.map(|j| bounds[j].0 + rng.next_f64() * (bounds[j].1 - bounds[j].0))
.collect()
})
.collect();
let mut fitness: Vec<f64> = pop.iter().map(|p| objective(p)).collect();
for t in 0..max_iter {
let best_idx = fitness
.iter()
.enumerate()
.min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
.unwrap()
.0;
let leader = pop[best_idx].clone();
let a = 2.0 * (1.0 - t as f64 / max_iter as f64);
for i in 0..pop_size {
let p = rng.next_f64();
if p < 0.5 {
let a_vec = 2.0 * a * rng.next_f64() - a;
if a_vec.abs() < 1.0 {
for j in 0..dim {
let d = (leader[j] - pop[i][j]).abs();
pop[i][j] = (leader[j] - a_vec * d).clamp(bounds[j].0, bounds[j].1);
}
} else {
let rand_idx = rng.next_usize(pop_size);
for (j, bj) in bounds.iter().enumerate() {
let d = (pop[rand_idx][j] - pop[i][j]).abs();
pop[i][j] = (pop[rand_idx][j] - a_vec * d).clamp(bj.0, bj.1);
}
}
} else {
let l = rng.next_f64() * 2.0 - 1.0;
let b_spiral = 1.0;
for j in 0..dim {
let d = (leader[j] - pop[i][j]).abs();
pop[i][j] =
(d * (b_spiral * l * std::f64::consts::PI).cos() * (b_spiral * l).exp()
+ leader[j])
.clamp(bounds[j].0, bounds[j].1);
}
}
fitness[i] = objective(&pop[i]);
}
}
let best = fitness
.iter()
.enumerate()
.min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
.unwrap()
.0;
pop[best].clone()
}