diff --git a/src/algorithms/ant_colony_tsp.rs b/src/algorithms/ant_colony_tsp.rs index ccadabd..e25f4b1 100644 --- a/src/algorithms/ant_colony_tsp.rs +++ b/src/algorithms/ant_colony_tsp.rs @@ -197,7 +197,7 @@ fn build_tour( for _ in 1..n { let current = *tour.last().unwrap(); // Build a probability vector over the unvisited candidates. - let mut probs: Vec<(usize, f64)> = (0..n) + let probs: Vec<(usize, f64)> = (0..n) .filter(|&j| !visited[j]) .map(|j| { let p = pheromone[current][j].max(0.0).powf(alpha) * eta[current][j].powf(beta); diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index ee50aac..3122a5f 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -17,6 +17,7 @@ pub mod random_search; pub mod simulated_annealing; pub mod spea2; pub mod tabu_search; +pub mod umda; pub use ant_colony_tsp::*; pub use cma_es::*; @@ -34,3 +35,4 @@ pub use random_search::*; pub use simulated_annealing::*; pub use spea2::*; pub use tabu_search::*; +pub use umda::*; diff --git a/src/algorithms/mopso.rs b/src/algorithms/mopso.rs index 8076216..d7e105f 100644 --- a/src/algorithms/mopso.rs +++ b/src/algorithms/mopso.rs @@ -4,7 +4,6 @@ use rand::Rng as _; use rand::seq::IndexedRandom; use crate::algorithms::parallel_eval::evaluate_batch; -use crate::core::candidate::Candidate; use crate::core::population::Population; use crate::core::problem::Problem; use crate::core::result::OptimizationResult; diff --git a/src/algorithms/umda.rs b/src/algorithms/umda.rs new file mode 100644 index 0000000..cb49c73 --- /dev/null +++ b/src/algorithms/umda.rs @@ -0,0 +1,289 @@ +//! `Umda` — Mühlenbein 1997 Univariate Marginal Distribution Algorithm for +//! binary (`Vec`) decisions. + +use rand::Rng as _; + +use crate::algorithms::parallel_eval::evaluate_batch; +use crate::core::candidate::Candidate; +use crate::core::objective::Direction; +use crate::core::population::Population; +use crate::core::problem::Problem; +use crate::core::result::OptimizationResult; +use crate::core::rng::rng_from_seed; +use crate::pareto::front::best_candidate; +use crate::traits::Optimizer; + +/// Configuration for [`Umda`]. +#[derive(Debug, Clone)] +pub struct UmdaConfig { + /// Sample size per generation. + pub population_size: usize, + /// Number of top members to use for the marginal estimate. + pub selected_size: usize, + /// Number of generations. + pub generations: usize, + /// Number of bits in each decision. + pub bits: usize, + /// Seed for the deterministic RNG. + pub seed: u64, +} + +impl Default for UmdaConfig { + fn default() -> Self { + Self { + population_size: 100, + selected_size: 50, + generations: 50, + bits: 32, + seed: 42, + } + } +} + +/// Univariate Marginal Distribution Algorithm. +/// +/// `Vec` decisions only; single-objective only. Each generation +/// estimates per-bit marginal probabilities from the top `selected_size` +/// members and samples the next population from the resulting independent +/// Bernoulli vector. Probabilities are clamped to +/// `[1 / (2·selected_size), 1 - 1 / (2·selected_size)]` (Laplace-style +/// smoothing) so the population never collapses to a single deterministic +/// string. +#[derive(Debug, Clone)] +pub struct Umda { + /// Algorithm configuration. + pub config: UmdaConfig, +} + +impl Umda { + /// Construct a `Umda`. + pub fn new(config: UmdaConfig) -> Self { + Self { config } + } +} + +impl

Optimizer

for Umda +where + P: Problem> + Sync, +{ + fn run(&mut self, problem: &P) -> OptimizationResult { + assert!(self.config.population_size >= 2, "Umda population_size must be >= 2"); + assert!( + self.config.selected_size >= 1, + "Umda selected_size must be >= 1", + ); + assert!( + self.config.selected_size <= self.config.population_size, + "Umda selected_size must be <= population_size", + ); + assert!(self.config.bits >= 1, "Umda bits must be >= 1"); + let objectives = problem.objectives(); + assert!( + objectives.is_single_objective(), + "Umda requires exactly one objective", + ); + let direction = objectives.objectives[0].direction; + let n = self.config.population_size; + let bits = self.config.bits; + let mu = self.config.selected_size; + let mut rng = rng_from_seed(self.config.seed); + + // Initial sample: uniform Bernoulli(0.5) across all bits. + let mut decisions: Vec> = (0..n) + .map(|_| (0..bits).map(|_| rng.random_bool(0.5)).collect()) + .collect(); + let mut population = evaluate_batch(problem, decisions.clone()); + let mut evaluations = population.len(); + + let smoothing = 1.0 / (2.0 * mu as f64); + let prob_min = smoothing; + let prob_max = 1.0 - smoothing; + + let mut best_seen: Option>> = None; + for c in &population { + let beats = match &best_seen { + None => true, + Some(b) => better_than_so(&c.evaluation, &b.evaluation, direction), + }; + if beats { + best_seen = Some(c.clone()); + } + } + + for _ in 0..self.config.generations { + // --- Phase 1: select top μ members --- + let mut order: Vec = (0..population.len()).collect(); + order.sort_by(|&a, &b| { + compare_so(&population[a].evaluation, &population[b].evaluation, direction) + }); + let selected: Vec<&Candidate>> = + order.iter().take(mu).map(|&i| &population[i]).collect(); + + // --- Phase 2: estimate per-bit marginals --- + let mut probs = vec![0.0_f64; bits]; + for c in &selected { + for (i, b) in c.decision.iter().enumerate() { + if *b { + probs[i] += 1.0; + } + } + } + for p in probs.iter_mut() { + *p = (*p / mu as f64).clamp(prob_min, prob_max); + } + + // --- Phase 3: sample a new population (uses RNG serially) --- + decisions = (0..n) + .map(|_| probs.iter().map(|&p| rng.random_bool(p)).collect()) + .collect(); + + // --- Phase 4: evaluate (parallel-friendly) --- + population = evaluate_batch(problem, decisions.clone()); + evaluations += population.len(); + + // Track best. + for c in &population { + let beats = match &best_seen { + None => true, + Some(b) => better_than_so(&c.evaluation, &b.evaluation, direction), + }; + if beats { + best_seen = Some(c.clone()); + } + } + } + + let best = best_seen.expect("at least one generation evaluated"); + let final_pop = vec![best.clone()]; + let front = vec![best.clone()]; + let best_opt = best_candidate(&final_pop, &objectives); + OptimizationResult::new( + Population::new(final_pop), + front, + best_opt, + evaluations, + self.config.generations, + ) + } +} + +fn compare_so( + a: &crate::core::evaluation::Evaluation, + b: &crate::core::evaluation::Evaluation, + direction: Direction, +) -> std::cmp::Ordering { + match (a.is_feasible(), b.is_feasible()) { + (true, false) => std::cmp::Ordering::Less, + (false, true) => std::cmp::Ordering::Greater, + (false, false) => a + .constraint_violation + .partial_cmp(&b.constraint_violation) + .unwrap_or(std::cmp::Ordering::Equal), + (true, true) => match direction { + Direction::Minimize => a.objectives[0] + .partial_cmp(&b.objectives[0]) + .unwrap_or(std::cmp::Ordering::Equal), + Direction::Maximize => b.objectives[0] + .partial_cmp(&a.objectives[0]) + .unwrap_or(std::cmp::Ordering::Equal), + }, + } +} + +fn better_than_so( + a: &crate::core::evaluation::Evaluation, + b: &crate::core::evaluation::Evaluation, + direction: Direction, +) -> bool { + compare_so(a, b, direction) == std::cmp::Ordering::Less +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::core::evaluation::Evaluation; + use crate::core::objective::{Objective, ObjectiveSpace}; + + /// OneMax: maximize the sum of true bits. + struct OneMax { + #[allow(dead_code)] + bits: usize, + } + impl Problem for OneMax { + type Decision = Vec; + + fn objectives(&self) -> ObjectiveSpace { + ObjectiveSpace::new(vec![Objective::maximize("bits")]) + } + + fn evaluate(&self, x: &Vec) -> Evaluation { + let count = x.iter().filter(|b| **b).count(); + Evaluation::new(vec![count as f64]) + } + } + + /// Trivial multi-objective problem to exercise the panic. + struct DummyMo; + impl Problem for DummyMo { + type Decision = Vec; + + fn objectives(&self) -> ObjectiveSpace { + ObjectiveSpace::new(vec![ + Objective::minimize("a"), + Objective::minimize("b"), + ]) + } + + fn evaluate(&self, _x: &Vec) -> Evaluation { + Evaluation::new(vec![0.0, 0.0]) + } + } + + #[test] + fn solves_onemax_20() { + let problem = OneMax { bits: 20 }; + let mut opt = Umda::new(UmdaConfig { + population_size: 50, + selected_size: 20, + generations: 30, + bits: 20, + seed: 1, + }); + let r = opt.run(&problem); + let best = r.best.unwrap(); + assert_eq!(best.evaluation.objectives[0], 20.0); + } + + #[test] + fn deterministic_with_same_seed() { + let problem = OneMax { bits: 16 }; + let cfg = UmdaConfig { + population_size: 30, + selected_size: 10, + generations: 10, + bits: 16, + seed: 99, + }; + let mut a = Umda::new(cfg.clone()); + let mut b = Umda::new(cfg); + let ra = a.run(&problem); + let rb = b.run(&problem); + assert_eq!( + ra.best.unwrap().evaluation.objectives, + rb.best.unwrap().evaluation.objectives, + ); + } + + #[test] + #[should_panic(expected = "exactly one objective")] + fn multi_objective_panics() { + let mut opt = Umda::new(UmdaConfig { + population_size: 10, + selected_size: 5, + generations: 1, + bits: 4, + seed: 0, + }); + let _ = opt.run(&DummyMo); + } +} diff --git a/src/prelude.rs b/src/prelude.rs index 23bdddc..7924eaa 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -27,5 +27,6 @@ pub use crate::algorithms::{ GeneticAlgorithm, GeneticAlgorithmConfig, HillClimber, HillClimberConfig, Ibea, IbeaConfig, Moead, MoeadConfig, Mopso, MopsoConfig, Nsga2, Nsga2Config, Nsga3, Nsga3Config, Paes, PaesConfig, ParticleSwarm, ParticleSwarmConfig, RandomSearch, RandomSearchConfig, SimulatedAnnealing, - SimulatedAnnealingConfig, Spea2, Spea2Config, TabuSearch, TabuSearchConfig, + SimulatedAnnealingConfig, Spea2, Spea2Config, TabuSearch, TabuSearchConfig, Umda, + UmdaConfig, };