genetic_algorithms 3.0.0

//! MOEA/D — Decomposition-based Multi-Objective Evolutionary Algorithm.
//!
//! ## Description
//!
//! MOEA/D (Zhang & Li 2007) decomposes a multi-objective problem into N scalar
//! sub-problems using weight vectors. Each sub-problem maintains its own current
//! solution and a **neighbourhood** of similar weight vectors. Offspring are
//! generated by selecting two parents from the neighbourhood, applying crossover
//! and mutation, and then competing only against neighbouring sub-problems via
//! **Tchebycheff** or **PBI** scalarization. An ideal point `z*` tracks the best
//! value found for each objective.
//!
//! Weight vectors are either auto-generated via the Das-Dennis simplex lattice
//! ([`MoeaDConfiguration::with_weight_vectors_auto`](configuration::MoeaDConfiguration::with_weight_vectors_auto))
//! or user-supplied
//! ([`MoeaDConfiguration::with_weight_vectors`](configuration::MoeaDConfiguration::with_weight_vectors)).
//!
//! Per generation, MOEA/D:
//! 1. Iterates all N sub-problems in order.
//! 2. For each sub-problem i:
//!    a. Randomly pick two parents from `neighbourhood[i]`.
//!    b. Apply crossover + mutation to produce one offspring.
//!    c. Evaluate offspring objectives.
//!    d. Update ideal point `z*` (per-objective minimum).
//!    e. Walk `neighbourhood[i]` — if the offspring's scalarised value is better
//!    than the current solution for neighbour j, replace it (up to
//!    `max_neighbor_replacements` replacements).
//! 3. Post-hoc non-dominated sort on the final population → ParetoFront.
//!
//! MOEA/D does NOT use Pareto ranking for selection — it is a decomposition-based
//! approach that is often faster than Pareto-based methods.
//!
//! ## When to Use
//!
//! - **Problem type:** Multi-objective (2+ objectives)
//! - **Variable type:** Continuous (real-valued), binary
//! - **Population structure:** Single population (one solution per sub-problem)
//! - **Key strength:** Computationally efficient — no expensive non-dominated
//!   sort inside the inner sub-problem loop (but a final post-hoc sort is run
//!   for the Pareto front output). Parallelisable via sub-problem distribution.
//! - **Key weakness:** Quality depends on the distribution of weight vectors.
//!   If the Pareto front does not align with the weight grid, coverage may be
//!   uneven.
//!
//! ## Quick Reference
//!
//! ### Mandatory Parameters
//!
//! | Parameter | Type | Default | Description |
//! |-----------|------|---------|-------------|
//! | `num_objectives` | `usize` | `3` | Number of objectives. |
//! | `population_size` | `usize` | `100` | Population size (one per sub-problem). |
//! | `max_generations` | `usize` | `200` | Maximum generations. |
//! | `init_fn` | `Fn` | — | Chromosome initialisation. |
//!
//! ### Optional Parameters
//!
//! | Parameter | Type | Default | Description |
//! |-----------|------|---------|-------------|
//! | `objective_directions` | `Vec<ObjectiveDirection>` | All `Minimize` | Per-objective Min/Max. |
//! | `scalarization` | `ScalarizationFn` | `Tchebycheff` | Tchebycheff or PBI with penalty θ. |
//! | `neighborhood_size` | `usize` | `20` | T — number of neighbours per sub-problem. |
//! | `max_neighbor_replacements` | `usize` | `2` | Max replacements per offspring. |
//! | `ga_config` | `GaConfiguration` | `Default` | GA operators, limits, RNG seed. |
//! | `observer` | `MoeaDObserver<U>` | `None` | Lifecycle observer. |
//!
//! ## Complete Example
//!
//! ```rust,no_run
//! // no_run: MOEA/D engine example — illustrative API usage, not a runnable benchmark
//! use genetic_algorithms::moead::MoeaDGa;
//! use genetic_algorithms::moead::configuration::{
//!     MoeaDConfiguration, ScalarizationFn,
//! };
//! use genetic_algorithms::configuration::GaConfiguration;
//!
//! let moead_config = MoeaDConfiguration::new()
//!     .with_num_objectives(3)
//!     .with_population_size(91)
//!     .with_max_generations(300)
//!     .with_weight_vectors_auto(12)
//!     .with_scalarization(ScalarizationFn::Tchebycheff)
//!     .with_neighborhood_size(20)
//!     .with_max_neighbor_replacements(2);
//!
//! let ga_config = GaConfiguration::default();
//! // let mut moead = MoeaDGa::<MyChromosome>::new(moead_config, ga_config)
//! //     .with_initialization_fn(|n, alleles, repeat| { /* ... */ })
//! //     .build()?;
//! //
//! // let pareto_front = moead.run()?;
//! // println!("Front size: {}", pareto_front.len());
//! ```
//!
//! ## Configuration Tips
//!
//! - The number of weight vectors (and thus sub-problems) is determined by
//!   `with_weight_vectors_auto(p)` or custom. The population_size should equal
//!   the number of weight vectors.
//! - `Tchebycheff` works well for most problems. Use `Pbi` with `theta = 5.0`
//!   when the Pareto front has sharp corners or is dominated by a single
//!   objective.
//! - `neighborhood_size` (T) controls exploitation vs exploration: small T →
//!   more local competition (faster convergence); large T → more information
//!   sharing (better diversity).
//! - The final Pareto front is extracted via post-hoc non-dominated sort — the
//!   algorithm itself runs decomposition-based selection.
//!
//! ## When to Choose This vs NSGA-III
//!
//! | Criterion | MOEA/D | NSGA-III |
//! |-----------|--------|----------|
//! | Mechanism | Weight-vector decomposition | Reference-point niche |
//! | Objectives | 2+ | 3+ (many-objective) |
//! | Selection | Scalarisation + neighbourhood | Niche preservation |
//! | Efficiency | No sort in sub-problem loop | Normalise + associate per gen |
//! | Pareto front | Even coverage if weights align | Uniform reference coverage |
//!
//! ## References
//!
//! - Zhang, Q., & Li, H. (2007). MOEA/D: A multiobjective evolutionary
//!   algorithm based on decomposition. _IEEE Trans. on Evolutionary
//!   Computation_, 11(6), 712–731.

pub mod configuration;

use crate::configuration::GaConfiguration;
use crate::error::GaError;
use crate::moead::configuration::{MoeaDConfiguration, ScalarizationFn};
use crate::multi_objective::non_dominated_sort::{
    assign_ranks, non_dominated_sort_with_directions,
};
use crate::multi_objective::pareto::{ParetoFront, ParetoIndividual};
use crate::observer::MoeaDObserver;
use crate::operations::{crossover, mutation};
use crate::traits::{InitializationFn, LinearChromosome, MutationOperator, VectorFitness};
use rand::Rng;
#[cfg(all(not(target_arch = "wasm32"), feature = "parallel"))]
use rayon::prelude::*;
use std::sync::Arc;
use std::time::Instant;

/// MOEA/D decomposition-based multi-objective genetic algorithm orchestrator.
///
/// # Type Parameters
///
/// * `U` - Chromosome type implementing `ChromosomeT`.
///
/// # Examples
///
/// ```rust,no_run
/// use genetic_algorithms::moead::MoeaDGa;
/// use genetic_algorithms::moead::configuration::MoeaDConfiguration;
/// use genetic_algorithms::configuration::GaConfiguration;
/// use genetic_algorithms::chromosomes::Range as RangeChromosome;
///
/// let moead_config = MoeaDConfiguration::default();
/// let ga_config = GaConfiguration::default();
/// let engine = MoeaDGa::<RangeChromosome<f64>>::new(moead_config, ga_config);
/// ```
pub struct MoeaDGa<U>
where
    U: LinearChromosome + VectorFitness,
{
    /// MOEA/D specific configuration.
    pub moead_config: MoeaDConfiguration,
    /// Base GA configuration (operators, limits).
    pub ga_config: GaConfiguration,
    /// Alleles template for initialization.
    pub alleles: Vec<U::Gene>,
    /// Initialization function.
    pub initialization_fn: Option<Arc<InitializationFn<U::Gene>>>,
    /// Optional structured lifecycle observer for MOEA/D-specific events.
    pub observer: Option<Arc<dyn MoeaDObserver<U> + Send + Sync>>,
}

impl<U> MoeaDGa<U>
where
    U: LinearChromosome + VectorFitness,
{
    /// Creates a new `MoeaDGa` with the given configurations.
    pub fn new(moead_config: MoeaDConfiguration, ga_config: GaConfiguration) -> Self {
        MoeaDGa {
            moead_config,
            ga_config,
            alleles: Vec::new(),
            initialization_fn: None,
            observer: None,
        }
    }

    /// Attaches a structured lifecycle observer that receives MOEA/D-specific hooks.
    pub fn with_observer(mut self, obs: Arc<dyn MoeaDObserver<U> + Send + Sync>) -> Self {
        self.observer = Some(obs);
        self
    }

    /// Dispatches an observer hook if an observer is attached. No-op when `self.observer` is `None`.
    #[inline]
    pub(crate) fn notify<F: FnOnce(&dyn MoeaDObserver<U>)>(&self, f: F) {
        if let Some(ref obs) = self.observer {
            f(obs.as_ref());
        }
    }

    /// Sets the alleles template.
    pub fn with_alleles(mut self, alleles: Vec<U::Gene>) -> Self {
        self.alleles = alleles;
        self
    }

    /// Sets the initialization function.
    pub fn with_initialization_fn<F>(mut self, f: F) -> Self
    where
        F: Fn(usize, Option<&[U::Gene]>) -> Vec<U::Gene> + Send + Sync + 'static,
    {
        self.initialization_fn = Some(Arc::new(f));
        self
    }

    /// Validates configuration and returns a ready-to-run instance.
    pub fn build(self) -> Result<Self, GaError> {
        self.validate()?;
        Ok(self)
    }

    /// Validates the MOEA/D configuration.
    ///
    /// # Errors
    ///
    /// Returns `GaError::InvalidMoeaDConfiguration` if parameters are invalid.
    pub fn validate(&self) -> Result<(), GaError> {
        if self.moead_config.num_objectives == 0 {
            return Err(GaError::InvalidMoeaDConfiguration(
                "num_objectives must be > 0".to_string(),
            ));
        }
        if self.moead_config.population_size < 2 {
            return Err(GaError::InvalidMoeaDConfiguration(
                "population_size must be >= 2".to_string(),
            ));
        }
        if self.initialization_fn.is_none() {
            return Err(GaError::InvalidMoeaDConfiguration(
                "initialization_fn is required".to_string(),
            ));
        }
        if !self.moead_config.objective_directions.is_empty()
            && self.moead_config.objective_directions.len() != self.moead_config.num_objectives
        {
            return Err(GaError::InvalidMoeaDConfiguration(format!(
                "objective_directions length ({}) must match num_objectives ({})",
                self.moead_config.objective_directions.len(),
                self.moead_config.num_objectives
            )));
        }
        // Das-Dennis subdivision count must be >= 1 to avoid a degenerate all-zero point.
        if let Some(p) = self.moead_config.weight_vectors_auto_p() {
            if p == 0 {
                return Err(GaError::InvalidMoeaDConfiguration(
                    "Das-Dennis subdivision count p must be >= 1".to_string(),
                ));
            }
        }
        // Weight vectors must be configured (auto or custom) per D-06.
        let wvs = self.moead_config.effective_weight_vectors();
        match wvs {
            None => {
                return Err(GaError::InvalidMoeaDConfiguration(
                    "weight vectors must be configured via with_weight_vectors_auto(p) or with_weight_vectors(vecs)".to_string(),
                ));
            }
            Some(vecs) => {
                if vecs.is_empty() {
                    return Err(GaError::InvalidMoeaDConfiguration(
                        "weight vector list must not be empty".to_string(),
                    ));
                }
                for (i, wv) in vecs.iter().enumerate() {
                    if wv.len() != self.moead_config.num_objectives {
                        return Err(GaError::InvalidMoeaDConfiguration(format!(
                            "weight vector {} has dimension {}, expected {}",
                            i,
                            wv.len(),
                            self.moead_config.num_objectives
                        )));
                    }
                }
            }
        }
        Ok(())
    }

    /// Validates configuration and returns the materialised weight vectors.
    ///
    /// Combines `validate()` and `effective_weight_vectors()` into a single
    /// call so `run()` does not invoke the Das-Dennis generator twice.
    pub(crate) fn validate_and_get_weight_vectors(&self) -> Result<Vec<Vec<f64>>, GaError> {
        if self.moead_config.num_objectives == 0 {
            return Err(GaError::InvalidMoeaDConfiguration(
                "num_objectives must be > 0".to_string(),
            ));
        }
        if self.moead_config.population_size < 2 {
            return Err(GaError::InvalidMoeaDConfiguration(
                "population_size must be >= 2".to_string(),
            ));
        }
        if self.initialization_fn.is_none() {
            return Err(GaError::InvalidMoeaDConfiguration(
                "initialization_fn is required".to_string(),
            ));
        }
        if !self.moead_config.objective_directions.is_empty()
            && self.moead_config.objective_directions.len() != self.moead_config.num_objectives
        {
            return Err(GaError::InvalidMoeaDConfiguration(format!(
                "objective_directions length ({}) must match num_objectives ({})",
                self.moead_config.objective_directions.len(),
                self.moead_config.num_objectives
            )));
        }
        if let Some(p) = self.moead_config.weight_vectors_auto_p() {
            if p == 0 {
                return Err(GaError::InvalidMoeaDConfiguration(
                    "Das-Dennis subdivision count p must be >= 1".to_string(),
                ));
            }
        }
        let vecs = self
            .moead_config
            .effective_weight_vectors()
            .ok_or_else(|| {
                GaError::InvalidMoeaDConfiguration(
                    "weight vectors must be configured via with_weight_vectors_auto(p) or with_weight_vectors(vecs)".to_string(),
                )
            })?;
        if vecs.is_empty() {
            return Err(GaError::InvalidMoeaDConfiguration(
                "weight vector list must not be empty".to_string(),
            ));
        }
        for (i, wv) in vecs.iter().enumerate() {
            if wv.len() != self.moead_config.num_objectives {
                return Err(GaError::InvalidMoeaDConfiguration(format!(
                    "weight vector {} has dimension {}, expected {}",
                    i,
                    wv.len(),
                    self.moead_config.num_objectives
                )));
            }
        }
        Ok(vecs)
    }
}

impl<U> MoeaDGa<U>
where
    U: LinearChromosome
        + VectorFitness
        + mutation::ValueMutable
        + crate::traits::RealValuedMutation,
{
    /// Runs the MOEA/D algorithm and returns the post-hoc Pareto front.
    ///
    /// Implements Zhang & Li 2007 Algorithm 1:
    /// 1. Validate and materialise weight vectors (single Das-Dennis call).
    /// 2. Precompute T nearest neighbours per sub-problem in weight-vector space.
    /// 3. Initialize population (one chromosome per sub-problem) and evaluate objectives.
    /// 4. Initialise ideal point z* from the starting population.
    /// 5. For each generation, iterate all N sub-problems. Per sub-problem i:
    ///    a. Sample two parents from neighbours\[i\] uniformly at random.
    ///    b. Apply crossover + mutation to produce one offspring chromosome.
    ///    c. Evaluate offspring objectives.
    ///    d. Update z* per component: z\*\[k\] = min(z\*\[k\], f_k(offspring)).
    ///    e. Walk neighbours\[i\] in order; for each j, if g(offspring, w_j, z*) <
    ///    g(pop\[j\].objectives, w_j, z*), replace population\[j\] with the offspring;
    /// 6. Post-hoc non-dominated sort over the final population; return rank-0 individuals as ParetoFront.
    ///
    /// # Errors
    ///
    /// - `GaError::InvalidMoeaDConfiguration` when configuration is incomplete.
    /// - `GaError::MutationError` when `Mutation::Differential` is configured (incompatible with MOEA/D's
    ///   per-sub-problem single-offspring loop, which lacks the population context Differential requires).
    /// - `GaError::CrossoverError` / `GaError::InitializationError` on operator failures.
    pub fn run(&mut self) -> Result<ParetoFront<U>, GaError> {
        // Single materialisation of weight vectors (avoids the double-Das-Dennis bug from NSGA-III WR-02).
        let weight_vectors = self.validate_and_get_weight_vectors()?;
        crate::rng::set_seed(self.ga_config.rng_seed);

        let pop_size = self.moead_config.population_size;
        let max_gens = self.moead_config.max_generations;
        let directions = self.moead_config.effective_directions();
        let scalarization = self.moead_config.scalarization;
        let num_objectives = self.moead_config.num_objectives;
        let max_replacements = self.moead_config.max_neighbor_replacements;

        // T capped at the population size to avoid out-of-bounds neighbour indices when
        // user supplies T > population_size (e.g., default T=20 with pop_size=15 in tests).
        let t_neigh = self
            .moead_config
            .neighborhood_size
            .min(weight_vectors.len())
            .max(1);

        // Step 2: precompute neighbourhoods.
        let neighbourhoods = precompute_neighbourhoods(&weight_vectors, t_neigh);

        // Step 3: initialise population (one individual per sub-problem; pop_size aligned with N).
        let mut population = self.initialize_population()?;

        // Runtime check: verify chromosome's fitness_values() matches num_objectives.
        if let Some(first) = population.first() {
            let got = first.chromosome.fitness_values().len();
            if got != self.moead_config.num_objectives {
                return Err(GaError::InvalidMoeaDConfiguration(format!(
                    "Expected {} objectives from fitness_values(), got {}",
                    self.moead_config.num_objectives, got
                )));
            }
        }

        // Step 4: initialise ideal point z* across all initial individuals.
        let mut ideal_point = vec![f64::INFINITY; num_objectives];
        for ind in &population {
            for (k, &f) in ind.objectives.iter().enumerate() {
                if f < ideal_point[k] {
                    ideal_point[k] = f;
                }
            }
        }
        // If any component is still INFINITY (empty population — defensive), fall back to 0.
        for v in ideal_point.iter_mut() {
            if !v.is_finite() {
                *v = 0.0;
            }
        }

        let n_subproblems = population.len();

        // Step 5: generation loop.
        for gen in 0..max_gens {
            let t_sort: Option<Instant> = if self.observer.is_some() {
                #[cfg(not(target_arch = "wasm32"))]
                {
                    Some(Instant::now())
                }
                #[cfg(target_arch = "wasm32")]
                {
                    None
                }
            } else {
                None
            };

            // Sub-problem update loop.
            let mut rng = crate::rng::make_rng();
            for i in 0..n_subproblems {
                // 5a: sample two parents from neighbours[i].
                let neigh = &neighbourhoods[i];
                let parent_a_idx = neigh[rng.random_range(0..neigh.len())];
                let parent_b_idx = neigh[rng.random_range(0..neigh.len())];

                // 5b: produce one offspring via crossover + mutation.
                let offspring_chrom = self.create_offspring_for_subproblem(
                    &population[parent_a_idx].chromosome,
                    &population[parent_b_idx].chromosome,
                    &mut rng,
                )?;

                // 5c: evaluate offspring objectives via VectorFitness.
                let mut offspring_chrom = offspring_chrom;
                offspring_chrom.calculate_fitness();
                let offspring_objectives = offspring_chrom.fitness_values().to_vec();

                // 5d: update ideal point z* incrementally per component.
                for (k, &f) in offspring_objectives.iter().enumerate() {
                    if f < ideal_point[k] {
                        ideal_point[k] = f;
                    }
                }

                // 5e: neighbourhood replacement, capped at max_replacements.
                let mut replacement_count = 0usize;
                for &j in neigh {
                    if replacement_count >= max_replacements {
                        break;
                    }
                    let g_offspring = scalarize(
                        &offspring_objectives,
                        &weight_vectors[i],
                        &ideal_point,
                        scalarization,
                    );
                    let g_current = scalarize(
                        &population[j].objectives,
                        &weight_vectors[j],
                        &ideal_point,
                        scalarization,
                    );
                    if g_offspring < g_current {
                        population[j] = ParetoIndividual::new(
                            offspring_chrom.clone(),
                            offspring_objectives.clone(),
                        );
                        replacement_count += 1;
                    }
                }
            }

            // Per-generation observer hooks: derive front_count from current population's ranks.
            // Compute fronts (used both for observer front_count and for stamping ranks).
            let obj_slices: Vec<&[f64]> = population
                .iter()
                .map(|ind| ind.objectives.as_slice())
                .collect();
            let fronts = non_dominated_sort_with_directions(&obj_slices, &directions);
            let mut ranks = vec![0usize; population.len()];
            assign_ranks(&mut ranks, &fronts);
            for (i, &r) in ranks.iter().enumerate() {
                population[i].rank = r;
            }
            let front_count = fronts.len();

            if let Some(start) = t_sort {
                self.notify(|obs| {
                    obs.on_non_dominated_sort_complete(gen, start.elapsed().as_secs_f64() * 1000.0)
                });
            }
            self.notify(|obs| obs.on_pareto_front_assigned(gen, front_count, population.len()));

            // Suppress pop_size dead-store warning when pop_size differs from n_subproblems.
            let _ = pop_size;
        }

        let front_individuals: Vec<ParetoIndividual<U>> =
            population.into_iter().filter(|ind| ind.rank == 0).collect();
        Ok(ParetoFront::new(front_individuals))
    }

    /// Initializes the population with random chromosomes and evaluates objectives in parallel.
    fn initialize_population(&self) -> Result<Vec<ParetoIndividual<U>>, GaError> {
        let init_fn = self.initialization_fn.as_ref().ok_or_else(|| {
            GaError::InitializationError("No initialization function set".to_string())
        })?;

        let pop_size = self.moead_config.population_size;
        let genes_per_chrom = match self.ga_config.limit_configuration.chromosome_length {
            crate::chromosomes::ChromosomeLength::Fixed(n) => n,
            crate::chromosomes::ChromosomeLength::Variable { .. } => {
                return Err(GaError::InvalidMoeaDConfiguration(
                    "ChromosomeLength::Variable is not yet supported (Phase 52). Use ChromosomeLength::Fixed.".into(),
                ));
            }
        };

        let alleles = if self.alleles.is_empty() {
            None
        } else {
            Some(self.alleles.as_slice())
        };

        let chromosomes: Vec<U> = crate::traits::initialize_chromosomes(
            pop_size,
            genes_per_chrom,
            alleles,
            init_fn,
            None,
            0,
        );

        #[cfg(all(not(target_arch = "wasm32"), feature = "parallel"))]
        let population: Vec<ParetoIndividual<U>> = chromosomes
            .into_par_iter()
            .map(|mut chrom| {
                chrom.calculate_fitness();
                let objectives = chrom.fitness_values().to_vec();
                ParetoIndividual::new(chrom, objectives)
            })
            .collect();
        #[cfg(any(target_arch = "wasm32", not(feature = "parallel")))]
        let population: Vec<ParetoIndividual<U>> = chromosomes
            .into_iter()
            .map(|mut chrom| {
                chrom.calculate_fitness();
                let objectives = chrom.fitness_values().to_vec();
                ParetoIndividual::new(chrom, objectives)
            })
            .collect();

        Ok(population)
    }

    /// Produces one offspring chromosome from two parents via crossover + mutation.
    ///
    /// Used only inside the MOEA/D sub-problem update loop. Differs from the
    /// NSGA-III batch `create_offspring()` because MOEA/D needs exactly one
    /// offspring per sub-problem per generation, not pop_size offspring per
    /// generation.
    fn create_offspring_for_subproblem(
        &self,
        parent_a: &U,
        parent_b: &U,
        rng: &mut impl Rng,
    ) -> Result<U, GaError> {
        let crossover_config = self.ga_config.crossover_configuration;
        let mutation_config = self.ga_config.mutation_configuration;
        let crossover_prob = crossover_config.probability_max.unwrap_or(1.0);
        let mut_prob = mutation_config.probability_max.unwrap_or(0.1);

        let p: f64 = rng.random();
        let mut children = if p <= crossover_prob {
            crossover::factory(parent_a, parent_b, crossover_config)?
        } else {
            vec![parent_a.clone(), parent_b.clone()]
        };

        // Pick a single offspring from the children vec (use first; defensive empty handling).
        let mut child = children.pop().unwrap_or_else(|| parent_a.clone());
        // Drop the rest deterministically (no use for the second child in MOEA/D).
        drop(children);

        // Mutation dispatch — trait-based single call (params are in the variant).
        let mp: f64 = rng.random();
        if mp <= mut_prob {
            if matches!(
                mutation_config.method,
                crate::operations::Mutation::Differential(..)
            ) {
                return Err(GaError::MutationError(
                    "Differential mutation is not supported in MOEA/D; \
                     use Cauchy, LevyFlight, Polynomial, or a standard mutation method instead."
                        .to_string(),
                ));
            }
            mutation_config
                .method
                .mutate(&mut child, &mutation_config.method)?;
        }

        Ok(child)
    }
}

/// Precomputes T nearest neighbours per sub-problem using Euclidean distance in weight-vector space.
///
/// `t` is the requested neighbourhood size; capped internally at `weight_vectors.len()`.
/// Each result vector starts with `i` itself (distance 0.0) and contains the t closest indices.
fn precompute_neighbourhoods(weight_vectors: &[Vec<f64>], t: usize) -> Vec<Vec<usize>> {
    let n = weight_vectors.len();
    let t = t.min(n).max(1);
    let mut neighbourhoods: Vec<Vec<usize>> = Vec::with_capacity(n);
    for i in 0..n {
        let mut dists: Vec<(usize, f64)> = (0..n)
            .map(|j| {
                let d_sq: f64 = weight_vectors[i]
                    .iter()
                    .zip(weight_vectors[j].iter())
                    .map(|(a, b)| (a - b).powi(2))
                    .sum();
                (j, d_sq.sqrt())
            })
            .collect();
        dists.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
        neighbourhoods.push(dists.into_iter().take(t).map(|(j, _)| j).collect());
    }
    neighbourhoods
}

/// Dispatches to Tchebycheff or PBI scalarization. See module docs for formulas.
fn scalarize(
    objectives: &[f64],
    weights: &[f64],
    ideal: &[f64],
    scalarization: ScalarizationFn,
) -> f64 {
    match scalarization {
        ScalarizationFn::Tchebycheff => objectives
            .iter()
            .zip(weights.iter())
            .zip(ideal.iter())
            .map(|((f_i, w_i), z_i)| w_i * (f_i - z_i).abs())
            .fold(f64::NEG_INFINITY, f64::max),
        ScalarizationFn::Pbi { theta } => {
            // d1 = signed projection of (f - z*) onto the unit weight direction.
            // d2 = perpendicular distance from (f - z*) to the line spanned by w.
            let w_norm_sq: f64 = weights.iter().map(|w| w * w).sum::<f64>().max(f64::EPSILON);
            let w_norm = w_norm_sq.sqrt();
            let d1: f64 = objectives
                .iter()
                .zip(ideal.iter())
                .zip(weights.iter())
                .map(|((f_i, z_i), w_i)| (f_i - z_i) * w_i)
                .sum::<f64>()
                / w_norm;
            let d2_sq: f64 = objectives
                .iter()
                .zip(ideal.iter())
                .zip(weights.iter())
                .map(|((f_i, z_i), w_i)| {
                    let diff = f_i - z_i;
                    let proj = d1 * w_i / w_norm;
                    (diff - proj).powi(2)
                })
                .sum();
            d1.abs() + theta * d2_sq.sqrt()
        }
    }
}