//! Multi-objective optimization with Pareto frontier analysis
//!
//! This module provides comprehensive multi-objective optimization including:
//! - Pareto frontier computation and maintenance
//! - Multi-objective optimization algorithms (NSGA-II, MOEA/D)
//! - Hypervolume and other quality indicators
//! - Trade-off analysis and visualization
//! - Interactive optimization with user preferences

use crate::{AutoMLError, AutoMLResult};
use rand::prelude::*;
use serde::{Deserialize, Serialize};
use std::cmp::Ordering;
use std::collections::HashMap;

/// Multi-objective solution with multiple objective values
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct MultiObjectiveSolution {
    pub id: String,
    pub parameters: HashMap<String, String>,
    pub objectives: Vec<f64>,
    pub objective_names: Vec<String>,
    pub constraints: Vec<f64>,
    pub metadata: HashMap<String, String>,
    pub dominance_rank: Option<usize>,
    pub crowding_distance: Option<f64>,
}

impl MultiObjectiveSolution {
    /// Create new multi-objective solution
    pub fn new(
        id: String,
        parameters: HashMap<String, String>,
        objectives: Vec<f64>,
        objective_names: Vec<String>,
    ) -> Self {
        Self {
            id,
            parameters,
            objectives,
            objective_names,
            constraints: Vec::new(),
            metadata: HashMap::new(),
            dominance_rank: None,
            crowding_distance: None,
        }
    }

    /// Check if this solution dominates another solution
    pub fn dominates(&self, other: &Self) -> bool {
        if self.objectives.len() != other.objectives.len() {
            return false;
        }

        let mut at_least_one_better = false;
        for i in 0..self.objectives.len() {
            if self.objectives[i] < other.objectives[i] {
                return false; // This solution is worse in at least one objective
            } else if self.objectives[i] > other.objectives[i] {
                at_least_one_better = true;
            }
        }

        at_least_one_better
    }

    /// Check if this solution is feasible (satisfies all constraints)
    pub fn is_feasible(&self) -> bool {
        self.constraints.iter().all(|&c| c <= 0.0) // Constraints are <= 0
    }

    /// Get objective value by name
    pub fn get_objective(&self, name: &str) -> Option<f64> {
        self.objective_names
            .iter()
            .position(|n| n == name)
            .map(|i| self.objectives[i])
    }

    /// Calculate distance to reference point
    pub fn distance_to_reference(&self, reference: &[f64]) -> f64 {
        if reference.len() != self.objectives.len() {
            return f64::INFINITY;
        }

        self.objectives
            .iter()
            .zip(reference)
            .map(|(obj, &ref_val)| (obj - ref_val).powi(2))
            .sum::<f64>()
            .sqrt()
    }
}

/// Pareto frontier with efficient operations
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct ParetoFrontier {
    solutions: Vec<MultiObjectiveSolution>,
    objective_names: Vec<String>,
    n_objectives: usize,
    reference_point: Option<Vec<f64>>,
    ideal_point: Option<Vec<f64>>,
    nadir_point: Option<Vec<f64>>,
}

impl ParetoFrontier {
    /// Create new empty Pareto frontier
    pub fn new(objective_names: Vec<String>) -> AutoMLResult<Self> {
        let n_objectives = objective_names.len();
        if n_objectives == 0 {
            return Err(AutoMLError::ConfigurationError(
                "Must have at least one objective".to_string(),
            ));
        }

        Ok(Self {
            solutions: Vec::new(),
            objective_names,
            n_objectives,
            reference_point: None,
            ideal_point: None,
            nadir_point: None,
        })
    }

    /// Add solution to frontier, updating Pareto set
    pub fn add_solution(&mut self, mut solution: MultiObjectiveSolution) -> AutoMLResult<bool> {
        if solution.objectives.len() != self.n_objectives {
            return Err(AutoMLError::ValidationError(
                "Solution has wrong number of objectives".to_string(),
            ));
        }

        solution.objective_names = self.objective_names.clone();

        // Check if solution is dominated by existing solutions
        for existing in &self.solutions {
            if existing.dominates(&solution) {
                return Ok(false); // Solution is dominated, not added
            }
        }

        // Remove solutions dominated by the new solution
        self.solutions
            .retain(|existing| !solution.dominates(existing));

        // Add the new solution
        self.solutions.push(solution);

        // Update ideal and nadir points
        self.update_reference_points();

        Ok(true)
    }

    /// Get all solutions on the frontier
    pub fn get_solutions(&self) -> &[MultiObjectiveSolution] {
        &self.solutions
    }

    /// Get number of solutions on frontier
    pub fn size(&self) -> usize {
        self.solutions.len()
    }

    /// Check if frontier is empty
    pub fn is_empty(&self) -> bool {
        self.solutions.is_empty()
    }

    /// Get best solution for specific objective
    pub fn get_best_for_objective(&self, objective_idx: usize) -> Option<&MultiObjectiveSolution> {
        if objective_idx >= self.n_objectives {
            return None;
        }

        self.solutions.iter().max_by(|a, b| {
            a.objectives[objective_idx]
                .partial_cmp(&b.objectives[objective_idx])
                .unwrap_or(Ordering::Equal)
        })
    }

    /// Get solution closest to reference point
    pub fn get_closest_to_reference(&self, reference: &[f64]) -> Option<&MultiObjectiveSolution> {
        if reference.len() != self.n_objectives {
            return None;
        }

        self.solutions.iter().min_by(|a, b| {
            a.distance_to_reference(reference)
                .partial_cmp(&b.distance_to_reference(reference))
                .unwrap_or(Ordering::Equal)
        })
    }

    /// Calculate hypervolume indicator
    pub fn calculate_hypervolume(&self, reference_point: Option<&[f64]>) -> AutoMLResult<f64> {
        if self.solutions.is_empty() {
            return Ok(0.0);
        }

        let ref_point = reference_point
            .map(<[f64]>::to_vec)
            .or_else(|| self.reference_point.clone())
            .unwrap_or_else(|| vec![0.0; self.n_objectives]);

        if ref_point.len() != self.n_objectives {
            return Err(AutoMLError::ValidationError(
                "Reference point has wrong dimensions".to_string(),
            ));
        }

        // Simplified hypervolume calculation for 2D case
        if self.n_objectives == 2 {
            self.calculate_hypervolume_2d(&ref_point)
        } else {
            // For higher dimensions, use Monte Carlo approximation
            self.calculate_hypervolume_monte_carlo(&ref_point, 100000)
        }
    }

    /// Calculate 2D hypervolume exactly
    fn calculate_hypervolume_2d(&self, reference: &[f64]) -> AutoMLResult<f64> {
        if self.solutions.is_empty() {
            return Ok(0.0);
        }

        // Sort solutions by first objective (descending)
        let mut sorted_solutions: Vec<_> = self.solutions.iter().collect();
        sorted_solutions.sort_by(|a, b| b.objectives[0].total_cmp(&a.objectives[0]));

        let mut hypervolume = 0.0;
        let mut prev_y = reference[1];

        for solution in sorted_solutions {
            if solution.objectives[1] > prev_y {
                let width = solution.objectives[0] - reference[0];
                let height = solution.objectives[1] - prev_y;

                if width > 0.0 && height > 0.0 {
                    hypervolume += width * height;
                    prev_y = solution.objectives[1];
                }
            }
        }

        Ok(hypervolume)
    }

    /// Calculate hypervolume using Monte Carlo sampling
    fn calculate_hypervolume_monte_carlo(
        &self,
        reference: &[f64],
        n_samples: usize,
    ) -> AutoMLResult<f64> {
        if self.solutions.is_empty() {
            return Ok(0.0);
        }

        // Find bounds for sampling
        let mut max_objectives = reference.to_vec();
        for solution in &self.solutions {
            for (i, &obj_val) in solution.objectives.iter().enumerate() {
                max_objectives[i] = max_objectives[i].max(obj_val);
            }
        }

        let mut rng = rand::thread_rng();
        let mut dominated_count = 0;

        for _ in 0..n_samples {
            // Generate random point in objective space
            let mut random_point = Vec::with_capacity(self.n_objectives);
            for i in 0..self.n_objectives {
                random_point.push(rng.gen_range(reference[i]..=max_objectives[i]));
            }

            // Check if point is dominated by any solution
            for solution in &self.solutions {
                let mut dominates = true;
                for j in 0..self.n_objectives {
                    if solution.objectives[j] < random_point[j] {
                        dominates = false;
                        break;
                    }
                }

                if dominates {
                    dominated_count += 1;
                    break;
                }
            }
        }

        // Calculate volume
        let mut total_volume = 1.0;
        for i in 0..self.n_objectives {
            total_volume *= max_objectives[i] - reference[i];
        }

        Ok(total_volume * (dominated_count as f64) / (n_samples as f64))
    }

    /// Calculate spacing metric (diversity measure)
    pub fn calculate_spacing(&self) -> f64 {
        if self.solutions.len() < 2 {
            return 0.0;
        }

        let mut distances = Vec::new();

        // Calculate minimum distance to other solutions for each solution
        for (i, sol1) in self.solutions.iter().enumerate() {
            let mut min_dist = f64::INFINITY;

            for (j, sol2) in self.solutions.iter().enumerate() {
                if i != j {
                    let dist = self.euclidean_distance(&sol1.objectives, &sol2.objectives);
                    min_dist = min_dist.min(dist);
                }
            }

            distances.push(min_dist);
        }

        // Calculate mean distance
        let mean_dist: f64 = distances.iter().sum::<f64>() / distances.len() as f64;

        // Calculate spacing (standard deviation of distances)
        let variance: f64 = distances
            .iter()
            .map(|&d| (d - mean_dist).powi(2))
            .sum::<f64>()
            / distances.len() as f64;

        variance.sqrt()
    }

    /// Calculate spread metric (extent of frontier)
    pub fn calculate_spread(&self) -> f64 {
        if self.solutions.is_empty() {
            return 0.0;
        }

        let mut total_spread = 0.0;

        for obj_idx in 0..self.n_objectives {
            let obj_values: Vec<f64> = self
                .solutions
                .iter()
                .map(|s| s.objectives[obj_idx])
                .collect();

            let min_val = obj_values.iter().fold(f64::INFINITY, |a, &b| a.min(b));
            let max_val = obj_values.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));

            total_spread += (max_val - min_val).powi(2);
        }

        total_spread.sqrt()
    }

    /// Merge with another Pareto frontier
    pub fn merge(&mut self, other: &Self) -> AutoMLResult<()> {
        if other.n_objectives != self.n_objectives {
            return Err(AutoMLError::ValidationError(
                "Frontiers have different number of objectives".to_string(),
            ));
        }

        for solution in &other.solutions {
            self.add_solution(solution.clone())?;
        }

        Ok(())
    }

    /// Filter solutions based on constraints
    pub fn filter_feasible(&mut self) {
        self.solutions.retain(MultiObjectiveSolution::is_feasible);
        self.update_reference_points();
    }

    /// Get knee point (best compromise solution)
    pub fn get_knee_point(&self) -> Option<&MultiObjectiveSolution> {
        if self.solutions.len() < 3 {
            return self.solutions.first();
        }

        let ideal = self.ideal_point.as_ref()?;
        let nadir = self.nadir_point.as_ref()?;

        // Normalize objectives
        let mut normalized_solutions = Vec::new();
        for solution in &self.solutions {
            let mut normalized_obj = Vec::new();
            for i in 0..self.n_objectives {
                let norm_val = if nadir[i] != ideal[i] {
                    (solution.objectives[i] - ideal[i]) / (nadir[i] - ideal[i])
                } else {
                    0.0
                };
                normalized_obj.push(norm_val);
            }
            normalized_solutions.push(normalized_obj);
        }

        // Find solution with maximum distance from line connecting extreme solutions
        let mut max_distance = 0.0;
        let mut knee_idx = 0;

        for (i, norm_obj) in normalized_solutions.iter().enumerate() {
            // Calculate distance to the diagonal line in normalized space
            let sum: f64 = norm_obj.iter().sum();
            let mean = sum / self.n_objectives as f64;
            let variance: f64 = norm_obj.iter().map(|&x| (x - mean).powi(2)).sum();
            let distance = variance.sqrt();

            if distance > max_distance {
                max_distance = distance;
                knee_idx = i;
            }
        }

        self.solutions.get(knee_idx)
    }

    /// Update ideal and nadir points
    fn update_reference_points(&mut self) {
        if self.solutions.is_empty() {
            self.ideal_point = None;
            self.nadir_point = None;
            return;
        }

        let mut ideal = vec![f64::INFINITY; self.n_objectives];
        let mut nadir = vec![f64::NEG_INFINITY; self.n_objectives];

        for solution in &self.solutions {
            for i in 0..self.n_objectives {
                ideal[i] = ideal[i].min(solution.objectives[i]);
                nadir[i] = nadir[i].max(solution.objectives[i]);
            }
        }

        self.ideal_point = Some(ideal);
        self.nadir_point = Some(nadir);
    }

    /// Calculate Euclidean distance between two points
    fn euclidean_distance(&self, point1: &[f64], point2: &[f64]) -> f64 {
        point1
            .iter()
            .zip(point2)
            .map(|(a, b)| (a - b).powi(2))
            .sum::<f64>()
            .sqrt()
    }

    /// Get ideal point (minimum values for each objective)
    pub fn get_ideal_point(&self) -> Option<&Vec<f64>> {
        self.ideal_point.as_ref()
    }

    /// Get nadir point (maximum values for each objective)
    pub fn get_nadir_point(&self) -> Option<&Vec<f64>> {
        self.nadir_point.as_ref()
    }

    /// Set reference point for hypervolume calculation
    pub fn set_reference_point(&mut self, reference: Vec<f64>) -> AutoMLResult<()> {
        if reference.len() != self.n_objectives {
            return Err(AutoMLError::ValidationError(
                "Reference point has wrong dimensions".to_string(),
            ));
        }
        self.reference_point = Some(reference);
        Ok(())
    }

    /// Export frontier data for visualization
    pub fn export_data(&self) -> HashMap<String, Vec<f64>> {
        let mut data = HashMap::new();

        for (i, obj_name) in self.objective_names.iter().enumerate() {
            let values: Vec<f64> = self.solutions.iter().map(|s| s.objectives[i]).collect();
            data.insert(obj_name.clone(), values);
        }

        data
    }
}

/// Multi-objective optimization algorithm implementations
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum MOOptimizationAlgorithm {
    /// NSGA-II (Non-dominated Sorting Genetic Algorithm II)
    NSGA2 {
        population_size: usize,
        n_generations: usize,
        crossover_rate: f64,
        mutation_rate: f64,
    },
    /// MOEA/D (Multi-Objective Evolutionary Algorithm based on Decomposition)
    MOEAD {
        population_size: usize,
        n_generations: usize,
        n_neighbors: usize,
        weight_vectors: Vec<Vec<f64>>,
    },
    /// SPEA2 (Strength Pareto Evolutionary Algorithm 2)
    SPEA2 {
        population_size: usize,
        archive_size: usize,
        n_generations: usize,
    },
}

/// Multi-objective optimization result
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct MOOptimizationResult {
    pub pareto_frontier: ParetoFrontier,
    pub all_evaluated_solutions: Vec<MultiObjectiveSolution>,
    pub hypervolume_history: Vec<f64>,
    pub spacing_history: Vec<f64>,
    pub n_evaluations: usize,
    pub optimization_time_seconds: f64,
}

/// Main multi-objective optimizer
pub struct MultiObjectiveOptimizer {
    algorithm: MOOptimizationAlgorithm,
    objective_names: Vec<String>,
    constraints: Vec<String>,
    random_state: Option<u64>,
}

impl MultiObjectiveOptimizer {
    /// Create new multi-objective optimizer
    pub fn new(
        algorithm: MOOptimizationAlgorithm,
        objective_names: Vec<String>,
    ) -> AutoMLResult<Self> {
        if objective_names.is_empty() {
            return Err(AutoMLError::ConfigurationError(
                "Must specify at least one objective".to_string(),
            ));
        }

        Ok(Self {
            algorithm,
            objective_names,
            constraints: Vec::new(),
            random_state: Some(42),
        })
    }

    /// Add constraint to optimization
    pub fn add_constraint(&mut self, constraint_name: String) {
        self.constraints.push(constraint_name);
    }

    /// Run multi-objective optimization
    pub async fn optimize<F>(&self, mut evaluation_fn: F) -> AutoMLResult<MOOptimizationResult>
    where
        F: FnMut(&HashMap<String, String>) -> AutoMLResult<(Vec<f64>, Vec<f64>)>,
    {
        let _start_time = std::time::Instant::now();
        let mut rng = StdRng::seed_from_u64(self.random_state.unwrap_or(42));

        match &self.algorithm {
            MOOptimizationAlgorithm::NSGA2 {
                population_size,
                n_generations,
                crossover_rate,
                mutation_rate,
            } => {
                self.run_nsga2(
                    &mut evaluation_fn,
                    *population_size,
                    *n_generations,
                    *crossover_rate,
                    *mutation_rate,
                    &mut rng,
                )
                .await
            }
            _ => {
                // Default implementation (simplified NSGA-II)
                self.run_nsga2(&mut evaluation_fn, 50, 100, 0.9, 0.1, &mut rng)
                    .await
            }
        }
    }

    /// Run NSGA-II algorithm
    async fn run_nsga2<F>(
        &self,
        evaluation_fn: &mut F,
        population_size: usize,
        n_generations: usize,
        _crossover_rate: f64,
        _mutation_rate: f64,
        rng: &mut StdRng,
    ) -> AutoMLResult<MOOptimizationResult>
    where
        F: FnMut(&HashMap<String, String>) -> AutoMLResult<(Vec<f64>, Vec<f64>)>,
    {
        let start_time = std::time::Instant::now();
        let mut all_solutions = Vec::new();
        let mut hypervolume_history = Vec::new();
        let mut spacing_history = Vec::new();
        let mut n_evaluations = 0;

        // Initialize population
        let mut population = Vec::new();
        for i in 0..population_size {
            // Generate random parameters (simplified)
            let mut parameters = HashMap::new();
            parameters.insert("param1".to_string(), rng.gen_range(0.0..1.0).to_string());
            parameters.insert("param2".to_string(), rng.gen_range(0.0..1.0).to_string());

            // Evaluate solution
            let (objectives, constraints) = evaluation_fn(&parameters)?;
            n_evaluations += 1;

            let solution = MultiObjectiveSolution {
                id: format!("gen0_ind{i}"),
                parameters,
                objectives,
                objective_names: self.objective_names.clone(),
                constraints,
                metadata: HashMap::new(),
                dominance_rank: None,
                crowding_distance: None,
            };

            population.push(solution.clone());
            all_solutions.push(solution);
        }

        // Evolution loop
        for generation in 0..n_generations {
            // Non-dominated sorting
            let fronts = self.non_dominated_sort(&population);

            // Calculate crowding distance for each front
            let mut ranked_population = Vec::new();
            for (rank, front) in fronts.iter().enumerate() {
                let mut front_with_distance = front.clone();
                self.calculate_crowding_distance(&mut front_with_distance);

                for mut solution in front_with_distance {
                    solution.dominance_rank = Some(rank);
                    ranked_population.push(solution);
                }
            }

            // Create new population through selection, crossover, and mutation
            let mut new_population = Vec::new();

            while new_population.len() < population_size {
                // Tournament selection
                let parent1 = self.tournament_selection(&ranked_population, rng);
                let _parent2 = self.tournament_selection(&ranked_population, rng);

                // Simple crossover and mutation (simplified)
                let mut child = parent1.clone();
                child.id = format!("gen{}_ind{}", generation + 1, new_population.len());

                // Mutate parameters slightly
                for value in child.parameters.values_mut() {
                    if rng.r#gen::<f64>() < 0.1 {
                        // 10% mutation rate
                        let current_val: f64 = value.parse().unwrap_or(0.5);
                        let mutated_val = (current_val + rng.gen_range(-0.1..0.1)).clamp(0.0, 1.0);
                        *value = mutated_val.to_string();
                    }
                }

                // Evaluate child
                let (objectives, constraints) = evaluation_fn(&child.parameters)?;
                n_evaluations += 1;

                child.objectives = objectives;
                child.constraints = constraints;
                child.dominance_rank = None;
                child.crowding_distance = None;

                new_population.push(child.clone());
                all_solutions.push(child);
            }

            population = new_population;

            // Calculate metrics for this generation
            let mut current_frontier = ParetoFrontier::new(self.objective_names.clone())?;
            for solution in &population {
                current_frontier.add_solution(solution.clone())?;
            }

            let hypervolume = current_frontier.calculate_hypervolume(None).unwrap_or(0.0);
            let spacing = current_frontier.calculate_spacing();

            hypervolume_history.push(hypervolume);
            spacing_history.push(spacing);
        }

        // Extract final Pareto frontier
        let mut final_frontier = ParetoFrontier::new(self.objective_names.clone())?;
        for solution in &all_solutions {
            final_frontier.add_solution(solution.clone())?;
        }

        let optimization_time = start_time.elapsed().as_secs_f64();

        Ok(MOOptimizationResult {
            pareto_frontier: final_frontier,
            all_evaluated_solutions: all_solutions,
            hypervolume_history,
            spacing_history,
            n_evaluations,
            optimization_time_seconds: optimization_time,
        })
    }

    /// Non-dominated sorting for NSGA-II
    fn non_dominated_sort(
        &self,
        population: &[MultiObjectiveSolution],
    ) -> Vec<Vec<MultiObjectiveSolution>> {
        let mut fronts = Vec::new();
        let mut domination_count = vec![0; population.len()];
        let mut dominated_solutions = vec![Vec::new(); population.len()];

        // Calculate domination relationships
        for i in 0..population.len() {
            for j in 0..population.len() {
                if i != j {
                    if population[i].dominates(&population[j]) {
                        dominated_solutions[i].push(j);
                    } else if population[j].dominates(&population[i]) {
                        domination_count[i] += 1;
                    }
                }
            }
        }

        // Find first front (non-dominated solutions)
        let mut current_front = Vec::new();
        for i in 0..population.len() {
            if domination_count[i] == 0 {
                current_front.push(population[i].clone());
            }
        }

        let mut front_index = 0;
        while !current_front.is_empty() {
            fronts.push(current_front.clone());

            let mut next_front = Vec::new();
            for sol_idx in 0..population.len() {
                if domination_count[sol_idx] == front_index + 1 {
                    // Check if this solution should be in the next front
                    let mut dominated_by_current_front = false;
                    for front_sol in &current_front {
                        if front_sol.dominates(&population[sol_idx]) {
                            dominated_by_current_front = true;
                            break;
                        }
                    }

                    if !dominated_by_current_front {
                        next_front.push(population[sol_idx].clone());
                    }
                }
            }

            current_front = next_front;
            front_index += 1;

            if front_index > population.len() {
                break; // Safety check
            }
        }

        fronts
    }

    /// Calculate crowding distance for solutions in a front
    fn calculate_crowding_distance(&self, front: &mut [MultiObjectiveSolution]) {
        if front.len() <= 2 {
            for solution in front.iter_mut() {
                solution.crowding_distance = Some(f64::INFINITY);
            }
            return;
        }

        // Initialize distances to 0
        for solution in front.iter_mut() {
            solution.crowding_distance = Some(0.0);
        }

        // Calculate distance for each objective
        for obj_idx in 0..self.objective_names.len() {
            // Sort by objective value
            front.sort_by(|a, b| {
                a.objectives[obj_idx]
                    .partial_cmp(&b.objectives[obj_idx])
                    .unwrap_or(Ordering::Equal)
            });

            // Set boundary solutions to infinite distance
            front[0].crowding_distance = Some(f64::INFINITY);
            front[front.len() - 1].crowding_distance = Some(f64::INFINITY);

            let obj_range =
                front[front.len() - 1].objectives[obj_idx] - front[0].objectives[obj_idx];

            if obj_range > 0.0 {
                for i in 1..front.len() - 1 {
                    let distance_increment = (front[i + 1].objectives[obj_idx]
                        - front[i - 1].objectives[obj_idx])
                        / obj_range;
                    let current_distance = front[i].crowding_distance.unwrap_or(0.0);
                    front[i].crowding_distance = Some(current_distance + distance_increment);
                }
            }
        }
    }

    /// Tournament selection for NSGA-II
    fn tournament_selection(
        &self,
        population: &[MultiObjectiveSolution],
        rng: &mut StdRng,
    ) -> MultiObjectiveSolution {
        let tournament_size = 2;
        let mut best = &population[rng.gen_range(0..population.len())];

        for _ in 1..tournament_size {
            let candidate = &population[rng.gen_range(0..population.len())];

            // Compare based on dominance rank and crowding distance
            if self.compare_solutions(candidate, best) {
                best = candidate;
            }
        }

        best.clone()
    }

    /// Compare two solutions for selection (rank first, then crowding distance)
    fn compare_solutions(&self, a: &MultiObjectiveSolution, b: &MultiObjectiveSolution) -> bool {
        let rank_a = a.dominance_rank.unwrap_or(usize::MAX);
        let rank_b = b.dominance_rank.unwrap_or(usize::MAX);

        if rank_a < rank_b {
            return true;
        } else if rank_a > rank_b {
            return false;
        }

        // Same rank, compare crowding distance (higher is better)
        let dist_a = a.crowding_distance.unwrap_or(0.0);
        let dist_b = b.crowding_distance.unwrap_or(0.0);

        dist_a > dist_b
    }
}