optirustic 1.2.3

use log::{debug, info, warn};
use rand::RngCore;
use std::collections::HashMap;
use std::fmt::{Display, Formatter};
use std::ops::Rem;
use std::path::PathBuf;

use crate::algorithms::nsga3::adaptive_ref_points::AdaptiveReferencePoints;
use crate::algorithms::nsga3::associate::AssociateToRefPoint;
use crate::algorithms::nsga3::niching::Niching;
use crate::algorithms::nsga3::normalise::Normalise;
use crate::algorithms::{Algorithm, NSGA2};
use crate::core::utils::get_rng;
use crate::core::{DataValue, Individual, OError};
use crate::operators::{
    Crossover, Mutation, ParetoConstrainedDominance, PolynomialMutation, PolynomialMutationArgs,
    Selector, SimulatedBinaryCrossover, SimulatedBinaryCrossoverArgs, TournamentSelector,
};
use crate::utils::{fast_non_dominated_sort, DasDarren1998, NumberOfPartitions};
use optirustic_macros::{as_algorithm, as_algorithm_args, impl_algorithm_trait_items};

#[cfg(feature = "python")]
use pyo3::exceptions::PyTypeError;
#[cfg(feature = "python")]
use pyo3::prelude::*;
#[cfg(feature = "python")]
use pyo3::IntoPyObjectExt;

mod adaptive_ref_points;
mod associate;
mod niching;
mod normalise;

/// The data key where the normalised objectives are stored for each [`Individual`].
const NORMALISED_OBJECTIVE_KEY: &str = "normalised_objectives";

/// The data key where the perpendicular distance to a reference point is stored for each [`Individual`].
const MIN_DISTANCE: &str = "distance";

/// The data key where the reference point with [`MIN_DISTANCE`] is stored for each [`Individual`].
const REF_POINT: &str = "reference_point";

/// The data key where the reference point index for [`REF_POINT`] is stored.
const REF_POINT_INDEX: &str = "reference_point_index";

/// The type for the number of individuals in the population.
#[derive(Serialize, Deserialize, Clone, Debug)]
pub enum Nsga3NumberOfIndividuals {
    /// The number of individuals are set equal to the number of reference points.
    EqualToReferencePointCount,
    /// Set a custom number of individuals. This must be larger than the number of reference points
    /// generated by setting the [`NSGA3Arg::number_of_partitions`].
    Custom(usize),
}

/// Convert a python object to `Nsga3NumberOfIndividuals`.
#[cfg(feature = "python")]
impl FromPyObject<'_> for Nsga3NumberOfIndividuals {
    fn extract_bound(obj: &Bound<'_, PyAny>) -> PyResult<Self> {
        if obj.is_none() {
            Ok(Nsga3NumberOfIndividuals::EqualToReferencePointCount)
        } else if let Ok(x) = obj.extract::<usize>() {
            Ok(Nsga3NumberOfIndividuals::Custom(x))
        } else {
            Err(PyTypeError::new_err("Invalid type".to_string()))
        }
    }
}

/// Convert `Nsga3NumberOfIndividuals` to PyObject to enable getter in `NSGA3Arg`.
#[cfg(feature = "python")]
impl<'py> IntoPyObject<'py> for Nsga3NumberOfIndividuals {
    type Target = PyAny;
    type Output = Bound<'py, Self::Target>;
    type Error = PyErr;

    fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
        let py_obj = match self {
            Nsga3NumberOfIndividuals::EqualToReferencePointCount => py.None().into_bound_py_any(py),
            Nsga3NumberOfIndividuals::Custom(total) => total.into_bound_py_any(py),
        };
        py_obj
    }
}

/// Input arguments for the NSGA3 algorithm.
#[as_algorithm_args]
#[cfg_attr(feature = "python", pyclass(get_all))]
pub struct NSGA3Arg {
    /// The number of individuals in the population.
    pub number_of_individuals: Nsga3NumberOfIndividuals,
    /// The number of partitions to use to calculate the reference points or weight.
    pub number_of_partitions: NumberOfPartitions,
    /// The options of the Simulated Binary Crossover (SBX) operator. This operator is used to
    /// generate new children by recombining the variables of parent solutions. This defaults to
    /// [`SimulatedBinaryCrossoverArgs::default()`].
    /// NOTE: it is advisable to use a large `distribution_index` to prevent the problem explained in
    /// Section IIa point #3 in the paper. With many objectives, "two distant parent solutions are
    /// likely to produce offspring solutions that are also distant from parents", which should be
    /// prevented.
    pub crossover_operator_options: Option<SimulatedBinaryCrossoverArgs>,
    /// The options to Polynomial Mutation (PM) operator used to mutate the variables of an
    /// individual. This defaults to [`PolynomialMutationArgs::default()`],
    /// with a distribution index or index parameter of `20` and variable probability equal `1`
    /// divided by the number of real variables in the problem (i.e., each variable will have the
    /// same probability of being mutated).
    pub mutation_operator_options: Option<PolynomialMutationArgs>,
    /// Instead of initialising the population with random variables, see the initial population
    /// with  the variable values from a JSON files exported with this tool. This option lets you
    /// restart the evolution from a previous generation; you can use any history file (exported
    /// when the field `export_history`) or the file exported when the stopping condition was reached.
    pub resume_from_file: Option<PathBuf>,
    /// The seed used in the random number generator (RNG). You can specify a seed in case you want
    /// to try to reproduce results. NSGA2 is a stochastic algorithm that relies on an RNG at
    /// different steps (when population is initially generated, during selection, crossover and
    /// mutation) and, as such, may lead to slightly different solutions. The seed is randomly
    /// picked if this is `None`.
    pub seed: Option<u64>,
}

/// Initialise the `NSGA3Arg` as python class.
#[cfg(feature = "python")]
#[pymethods]
impl NSGA3Arg {
    #[new]
    #[pyo3(signature = (number_of_individuals, number_of_partitions, stopping_condition, crossover_operator_options=None, mutation_operator_options=None, resume_from_file=None, parallel=None, export_history=None, seed=None))]
    fn py_new(
        number_of_individuals: PyObject,
        number_of_partitions: NumberOfPartitions,
        stopping_condition: StoppingCondition,
        crossover_operator_options: Option<SimulatedBinaryCrossoverArgs>,
        mutation_operator_options: Option<PolynomialMutationArgs>,
        resume_from_file: Option<PathBuf>,
        parallel: Option<bool>,
        export_history: Option<ExportHistory>,
        seed: Option<u64>,
    ) -> PyResult<Self> {
        let number_of_individuals = Python::with_gil(|py| number_of_individuals.extract(py))?;
        Ok(NSGA3Arg {
            number_of_individuals,
            number_of_partitions,
            crossover_operator_options,
            mutation_operator_options,
            seed,
            stopping_condition,
            resume_from_file,
            parallel,
            export_history,
        })
    }

    pub fn __repr__(&self) -> PyResult<String> {
        Ok(format!(
            "NSGA3Arg(number_of_individuals={:?}, number_of_partitions={:?}, stopping_condition={})",
            self.number_of_individuals, self.number_of_partitions, self.stopping_condition
        ))
    }

    fn __str__(&self) -> String {
        self.__repr__().unwrap()
    }
}

/// The Non-dominated Sorting Genetic Algorithm (NSGA3).
///
/// Implemented based on:
/// > K. Deb and H. Jain, "An Evolutionary Many-Objective Optimization Algorithm Using
/// > Reference-Point-Based Non-dominated Sorting Approach, Part I: Solving Problems With Box
/// > Constraints," in IEEE Transactions on Evolutionary Computation, vol. 18, no. 4, pp. 577-601,
/// > Aug. 2014, doi: 10.1109/TEVC.2013.2281535
///
/// See: <https://10.1109/TEVC.2013.2281535>
/// # Example - solve the DTLZ1 problem
/// ```rust
#[doc = include_str!("../../../examples/nsga3_dtlz1.rs")]
/// ```
#[as_algorithm(NSGA3Arg)]
pub struct NSGA3 {
    /// The vector of reference points
    reference_points: Vec<Vec<f64>>,
    /// The number of original reference points. This equals the size of `reference_points` if
    /// the `adaptive` option is set to `false`.
    number_of_or_reference_points: usize,
    /// The gap between the reference points.
    ref_point_gap: f64,
    /// The ideal point coordinates when the algorithm starts up to the current evolution
    ideal_point: Vec<f64>,
    /// The operator to use to select the individuals for reproduction. This is a binary tournament
    /// selector ([`TournamentSelector`]) with the [`ParetoConstrainedDominance`] comparison operator.
    selector_operator: TournamentSelector<ParetoConstrainedDominance>,
    /// The SBX operator to use to generate a new children by recombining the variables of parent
    /// solutions.
    crossover_operator: SimulatedBinaryCrossover,
    /// The PM operator to use to mutate the variables of an individual.
    mutation_operator: PolynomialMutation,
    /// The seed to use.
    rng: Box<dyn RngCore>,
    /// Whether to use the reference point adaptive approach and convert the algorithm from `NSGA3`
    /// to `A-NSGA3`.
    adaptive: bool,
}

impl NSGA3 {
    /// Initialise the NSGA3 algorithm.
    ///
    /// # Arguments
    ///
    /// * `problem`: The problem being solved.
    /// * `options`: The [`NSGA3Arg`] arguments to customise the algorithm behaviour.
    /// * `adaptive`: Whether to use the adaptive approach to handle the reference points. Normally
    /// reference points are fixed and never change position, and `NSGA3` will try to associate one
    /// individual to each reference point. However, there are some problems where the reference lines
    /// never intersects the Pareto front and no Pareto-optimal solutions are associated to those
    /// reference points. The reference points intersecting the front may have too many points
    /// associated and create crowded solutions.
    /// To improve the solutions for these kind of problems, when this option is set to `true`, new
    /// reference points will adaptively be added around the provided points to reduce crowding and
    /// force one solution to be associated to preferably only on reference direction.
    /// This option converts `NSGA3` to [`crate::algorithms::AdaptiveNSGA3`] (where `A` stands for
    /// adaptive) and it is explained in Section VII of Jain and Deb (2013). For a detailed
    /// explanation about the implementation see `AdaptiveReferencePoints`.
    ///
    /// returns: `NSGA3`.
    pub fn new(problem: Problem, options: NSGA3Arg, adaptive: bool) -> Result<Self, OError> {
        let name = if !adaptive {
            "NSGA3".to_string()
        } else {
            "AdaptiveNSGA3".to_string()
        };
        let nsga3_args = options.clone();

        let das_darren = DasDarren1998::new(
            problem.number_of_objectives(),
            &options.number_of_partitions,
        )?;
        let reference_points = das_darren.get_weights();
        let number_of_or_reference_points = reference_points.len();

        info!(
            "Created {} reference directions",
            number_of_or_reference_points
        );

        // create the population
        let mut number_of_individuals = match options.number_of_individuals {
            Nsga3NumberOfIndividuals::EqualToReferencePointCount => reference_points.len(),
            Nsga3NumberOfIndividuals::Custom(count) => {
                if count < 3 {
                    return Err(OError::AlgorithmInit(
                        name,
                        "The population size must have at least 3 individuals".to_string(),
                    ));
                }
                // add tolerance (for example Deb et al. sometimes uses pop + 2 = ref_point)
                if count - 2 > das_darren.number_of_points() as usize {
                    return Err(OError::AlgorithmInit(
                        name,
                        format!(
                            concat!(
                            "The number of individuals ({}) must be larger than the number of ",
                            "reference points ({}) to prevent unexpected behaviours. It is always ",
                            "suggested to associate at least one individual to a reference point"),
                            count,
                            das_darren.number_of_points()
                        ),
                    ));
                }
                count
            }
        };
        debug!("Population size set to {}", number_of_individuals);

        // force the population size as multiple of 2 so that the new number of generated offsprings
        // matches `number_of_individuals`
        if number_of_individuals.rem(2) != 0 {
            number_of_individuals -= 1;
            warn!(
                "The population size was reduced to {} so that it is a multiple of 2",
                number_of_individuals
            );
        }

        let problem = Arc::new(problem);
        let population = if let Some(init_file) = options.resume_from_file {
            info!("Loading initial population from {:?}", init_file);
            NSGA3::seed_population_from_file(
                problem.clone(),
                &name,
                number_of_individuals,
                &init_file,
            )?
        } else {
            info!("Created initial random population");
            Population::init(problem.clone(), number_of_individuals)
        };

        let selector_operator = TournamentSelector::<ParetoConstrainedDominance>::new(2);
        let mutation_options = match options.mutation_operator_options {
            Some(o) => o,
            None => PolynomialMutationArgs::default(problem.clone().as_ref()),
        };
        let mutation_operator = PolynomialMutation::new(mutation_options.clone())?;

        let crossover_options = options.crossover_operator_options.unwrap_or_default();
        let crossover_operator = SimulatedBinaryCrossover::new(crossover_options.clone())?;

        info!(
            "{}",
            NSGA2::algorithm_option_str(&problem, &crossover_options, &mutation_options)
        );

        Ok(Self {
            number_of_individuals,
            reference_points,
            number_of_or_reference_points,
            ref_point_gap: das_darren.gap(),
            ideal_point: vec![f64::INFINITY; problem.number_of_objectives()],
            population,
            problem,
            selector_operator,
            crossover_operator,
            mutation_operator,
            generation: 0,
            nfe: 0,
            stopping_condition: options.stopping_condition,
            start_time: Instant::now(),
            parallel: options.parallel.unwrap_or(true),
            export_history: options.export_history,
            rng: get_rng(options.seed),
            args: nsga3_args,
            adaptive,
        })
    }

    /// Get the normalised objective data stored in the `individual`.
    ///
    /// # Arguments
    ///
    /// * `individual`: The individual reference with the data.
    ///
    /// returns: `Result<DataValue, OError>`
    fn get_normalised_objectives(individual: &Individual) -> Result<DataValue, OError> {
        individual.get_data(NORMALISED_OBJECTIVE_KEY)
    }

    /// For each reference point, count selected individuals in P_{t+1} associated with it. This
    /// returns `rho_j` which is a lookup map, mapping the reference point index to the number of
    /// linked individuals.
    ///
    /// # Arguments
    ///
    /// * `selected_individuals`: The individuals to use to count the association to the reference
    ///    points.
    /// * `reference_points`: The reference points.
    ///
    /// returns: `Result<HashMap<usize, usize>, OError>`
    fn get_association_map(
        selected_individuals: &Population,
        reference_points: &[Vec<f64>],
    ) -> Result<HashMap<usize, usize>, OError> {
        let mut rho_j: HashMap<usize, usize> = HashMap::new();
        for ind in selected_individuals.individuals() {
            let ref_point_index = ind.get_data(REF_POINT_INDEX);
            match ref_point_index {
                Ok(index) => {
                    let index = index.as_usize()?;
                    rho_j.entry(index).and_modify(|v| *v += 1).or_insert(1);
                }
                Err(_) => continue,
            }
        }
        // fill the rest
        for ref_point_index in 0..reference_points.len() {
            rho_j.entry(ref_point_index).or_insert(0);
        }

        Ok(rho_j)
    }

    /// Get the reference points used in the evolution.
    ///
    /// return: `Vec<Vec<f64>>`
    pub fn reference_points(&self) -> Vec<Vec<f64>> {
        self.reference_points.clone()
    }
}

/// Implementation of Section IV of the paper.
#[impl_algorithm_trait_items(NSGA3Arg)]
impl Algorithm<NSGA3Arg> for NSGA3 {
    /// This assesses the initial random population.
    ///
    /// return: `Result<(), OError>`
    fn initialise(&mut self) -> Result<(), OError> {
        info!("Evaluating initial population");
        if self.parallel {
            NSGA3::do_parallel_evaluation(self.population.individuals_as_mut(), &mut self.nfe)?;
        } else {
            NSGA3::do_evaluation(self.population.individuals_as_mut(), &mut self.nfe)?;
        }

        info!("Initial evaluation completed");
        self.generation += 1;

        Ok(())
    }

    /// Evolve the population. The first part of this code comes from NSGA2::evolve(). NSGA3 mainly
    /// differs in the survival method.
    fn evolve(&mut self) -> Result<(), OError> {
        // Create the new population, based on the population at the previous time-step, of size
        // self.number_of_individuals. The loop adds two individuals at the time.
        debug!("Generating new population (selection + crossover + mutation)");
        let mut offsprings: Vec<Individual> = Vec::new();
        for _ in 0..self.number_of_individuals / 2 {
            let parents =
                self.selector_operator
                    .select(self.population.individuals(), 2, &mut self.rng)?;

            // generate the 2 children with crossover
            let children = self.crossover_operator.generate_offsprings(
                &parents[0],
                &parents[1],
                &mut self.rng,
            )?;

            // mutate them
            offsprings.push(
                self.mutation_operator
                    .mutate_offspring(&children.child1, &mut self.rng)?,
            );
            offsprings.push(
                self.mutation_operator
                    .mutate_offspring(&children.child2, &mut self.rng)?,
            );
        }
        debug!("Combining parents and offsprings in new population");
        self.population.add_new_individuals(offsprings);
        debug!("New population size is {}", self.population.len());

        debug!("Evaluating population");
        if self.parallel {
            NSGA3::do_parallel_evaluation(self.population.individuals_as_mut(), &mut self.nfe)?;
        } else {
            NSGA3::do_evaluation(self.population.individuals_as_mut(), &mut self.nfe)?;
        }
        debug!("Evaluation done");

        debug!("Calculating fronts and ranks for new population");
        let sorting_results = fast_non_dominated_sort(self.population.individuals_as_mut(), false)?;
        debug!("Collected {} fronts", sorting_results.fronts.len());

        debug!("Selecting best individuals");
        // this is S_t in the paper, the population with the last front
        let mut new_population = Population::new();

        // Algorithm 1 in paper, step 5-7 - Fill the new population up to the last front.
        let mut last_front: Option<Vec<Individual>> = None;
        for (fi, front) in sorting_results.fronts.into_iter().enumerate() {
            if new_population.len() + front.len() <= self.number_of_individuals {
                // population does not overflow with new front
                debug!("Adding front #{} (size: {})", fi + 1, front.len());
                new_population.add_new_individuals(front);
            } else if new_population.len() == self.number_of_individuals {
                debug!("Population reached target size");
                break;
            } else {
                // Algorithm 1, step 12. Population is filled up to front l-1 (P_{t+1})
                debug!(
                    "Population almost full ({} individuals)",
                    new_population.len()
                );
                // this is F_l
                last_front = Some(front);
                break;
            }
        }

        // Algorithm 1, step 11-18 - Complete the population using the members from the last front.
        if let Some(last_front) = last_front {
            // store the last population index containing individuals up to front F_l
            let first_dom_index = new_population.len();
            let missing_item_count = self.number_of_individuals - new_population.len();
            debug!("{missing_item_count} must be added from the last front");

            // add the last_front F_l to create S_t
            new_population.add_new_individuals(last_front);

            // Algorithm 1, step 14 - Calculate f_n
            debug!("Normalising all individuals");
            let mut norm =
                Normalise::new(&mut self.ideal_point, new_population.individuals_as_mut())?;
            norm.calculate()?;

            // Algorithm 1, step 15
            debug!("Associating reference points to all individuals");
            let mut assoc = AssociateToRefPoint::new(
                new_population.individuals_as_mut(),
                &self.reference_points,
            )?;
            assoc.calculate()?;

            // Algorithm 1, step 16
            // re-split population in P_{t+1} (S_t without the last front) and individuals in front F_l
            let mut potential_individuals = new_population.drain(first_dom_index..);
            // rename variable for clarity
            let mut selected_individuals = new_population;

            // for each reference point count selected individuals in P_{t+1} associated with it.
            // rho_j is a lookup map mapping the reference point index to the number of linked
            // individuals
            let mut rho_j =
                NSGA3::get_association_map(&selected_individuals, &self.reference_points)?;

            // Algorithm 4 - Niching
            debug!("Niching");
            let mut n = Niching::new(
                &mut selected_individuals,
                &mut potential_individuals,
                missing_item_count,
                &mut rho_j,
                &mut self.rng,
            )?;
            n.calculate()?;

            // update the population
            self.population = selected_individuals;

            // add new refernece points
            if self.adaptive {
                let mut a = AdaptiveReferencePoints::new(
                    &mut self.reference_points,
                    &mut rho_j,
                    self.number_of_or_reference_points,
                    self.ref_point_gap,
                )?;
                a.calculate()?;
            }
        } else {
            // update the population
            self.population = new_population;
        }

        self.generation += 1;
        Ok(())
    }

    fn additional_export_data(&self) -> Option<HashMap<String, DataValue>> {
        let mut data = HashMap::new();
        let mut points: Vec<DataValue> = Vec::new();
        for point in self.reference_points() {
            points.push(DataValue::Vector(point));
        }
        data.insert(
            "reference_points".to_string(),
            DataValue::DataVector(points),
        );
        data.insert(
            "ideal_point".to_string(),
            DataValue::Vector(self.ideal_point.clone()),
        );
        Some(data)
    }
}

#[cfg(test)]
mod test_problems {
    use float_cmp::{approx_eq, assert_approx_eq};

    use optirustic_macros::test_with_retries;

    use crate::algorithms::{
        Algorithm, NSGA3Arg, Nsga3NumberOfIndividuals, StoppingCondition, NSGA3,
    };
    use crate::core::builtin_problems::{DTLZ1Problem, DTLZ2Problem};
    use crate::core::test_utils::{assert_approx_array_eq, check_value_in_range};
    use crate::operators::{PolynomialMutationArgs, SimulatedBinaryCrossoverArgs};
    use crate::utils::{NumberOfPartitions, TwoLayerPartitions};

    /// Test the DTLZ1 problem from Deb et al. (2013)
    fn test_dtlz1(number_objectives: usize, max_gen: u32) {
        // see Table I
        let k: usize = 5;
        let number_variables: usize = number_objectives + k - 1; // M + k - 1 with k = 5 (Section Va)
        let problem = DTLZ1Problem::create(number_variables, number_objectives, false).unwrap();
        // The number of partitions used in the paper when from section 5
        let number_of_partitions = match number_objectives {
            3 => NumberOfPartitions::OneLayer(12),
            5 => NumberOfPartitions::OneLayer(6),
            8 | 10 => NumberOfPartitions::TwoLayers(TwoLayerPartitions {
                boundary_layer: 3,
                inner_layer: 2,
                scaling: None,
            }),
            15 => NumberOfPartitions::TwoLayers(TwoLayerPartitions {
                boundary_layer: 2,
                inner_layer: 1,
                scaling: None,
            }),
            _ => panic!("Objective count not supported"),
        };
        // number of individuals - from Table I
        let pop_size: usize = match number_objectives {
            3 => 92,
            5 => 212,
            8 => 156,
            10 => 276,
            15 => 136,

            _ => panic!("Objective count not supported"),
        };
        let expected_ref_points: usize = match number_objectives {
            3 => 91,
            5 => 210,
            8 => 156,
            10 => 275,
            15 => 135,

            _ => panic!("Objective count not supported"),
        };

        // see Table II
        let crossover_operator_options = SimulatedBinaryCrossoverArgs {
            distribution_index: 30.0,
            ..SimulatedBinaryCrossoverArgs::default()
        };
        // eta_m = 20 - probability  1/n_vars
        let mutation_operator_options = PolynomialMutationArgs::default(&problem);

        let args = NSGA3Arg {
            // see Table I
            number_of_individuals: Nsga3NumberOfIndividuals::Custom(pop_size),
            number_of_partitions,
            crossover_operator_options: Some(crossover_operator_options),
            mutation_operator_options: Some(mutation_operator_options),
            // see Table III
            stopping_condition: StoppingCondition::MaxGeneration(max_gen),
            parallel: None,
            resume_from_file: None,
            export_history: None,
            seed: Some(1),
        };

        let mut algo = NSGA3::new(problem, args, false).unwrap();
        assert_eq!(algo.reference_points().len(), expected_ref_points);

        algo.run().unwrap();
        let results = algo.get_results();

        let expected_vars = vec![0.5; number_variables];
        let mut invalid_individuals: usize = 0;
        for ind in results.individuals {
            // All objective points lie on the plane passing through the 0.5 intercept on each axis (i.e.
            // the sum of the objective coordinate is close to 0.5). Because of randomness a few solutions
            // may breach this condition.
            let obj_sum: f64 = ind.get_objective_values().unwrap().iter().sum();
            let outside_range_data = approx_eq!(f64, obj_sum, 0.5, epsilon = 0.01);
            if !outside_range_data {
                invalid_individuals += 1;
            }

            // All variables in x_M must be 0.5
            let vars: Vec<f64> = ((number_variables - k + 1)..=number_variables)
                .map(|i| {
                    ind.get_variable_value(format!("x{i}").as_str())
                        .unwrap()
                        .as_real()
                        .unwrap()
                })
                .collect();
            assert_approx_array_eq(&vars, &expected_vars, Some(0.01));
        }

        // about 90% of solutions are ideal
        if invalid_individuals > 10 {
            panic!("Found {invalid_individuals} individuals not meeting the ideal solution");
        }
    }

    #[test_with_retries(10)]
    /// Test the DTLZ1 problem with M=3 and MaxGeneration = 400 (Table III of NSGA3 paper)
    fn test_dtlz1_obj_3() {
        test_dtlz1(3, 400);
    }

    #[test_with_retries(5)]
    /// Test the DTLZ1 problem with M=5 and MaxGeneration = 600 (Table III of NSGA3 paper)
    fn test_dtlz1_obj_5() {
        test_dtlz1(5, 600);
    }

    #[test_with_retries(5)]
    /// Test the DTLZ1 problem with M=8 and MaxGeneration = 750 (Table III of NSGA3 paper)
    fn test_dtlz1_obj_8() {
        test_dtlz1(8, 750);
    }

    // These two tests take too long to run on pipeline
    // #[test_with_retries(3)]
    // /// Test the DTLZ1 problem with M=10 and MaxGeneration = 1000 (Table III of NSGA3 paper)
    // fn test_dtlz1_obj_10() {
    //     test_dtlz1(10, 1000);
    // }

    // #[test_with_retries(3)]
    // /// Test the DTLZ1 problem with M=15 and MaxGeneration = 1500 (Table III of NSGA3 paper)
    // fn test_dtlz1_obj_15() {
    //     test_dtlz1(15, 1500);
    // }

    #[test_with_retries(10)]
    /// Test the DTLZ2 problem from Deb et al. (2013) with M=3 (see Table III).
    fn test_dtlz2() {
        // see Table I
        let number_objectives: usize = 3;
        let k: usize = 10;
        let number_variables: usize = number_objectives + k - 1; // M + k - 1 with k = 5 (Section Va)
        let problem = DTLZ2Problem::create(number_variables, number_objectives).unwrap();
        // The number of partitions used in the paper when M=3 - Table I
        let number_of_partitions = NumberOfPartitions::OneLayer(12);

        // see Table II
        let crossover_operator_options = SimulatedBinaryCrossoverArgs {
            distribution_index: 30.0,
            ..SimulatedBinaryCrossoverArgs::default()
        };
        // eta_m = 20 - probability  1/n_vars
        let mutation_operator_options = PolynomialMutationArgs::default(&problem);

        let args = NSGA3Arg {
            // see Table I
            number_of_individuals: Nsga3NumberOfIndividuals::Custom(92),
            number_of_partitions,
            crossover_operator_options: Some(crossover_operator_options),
            mutation_operator_options: Some(mutation_operator_options),
            // see Table III
            stopping_condition: StoppingCondition::MaxGeneration(400),
            resume_from_file: None,
            parallel: None,
            export_history: None,
            seed: Some(1),
        };

        let mut algo = NSGA3::new(problem, args, false).unwrap();
        assert_eq!(algo.reference_points().len(), 91);

        algo.run().unwrap();
        let results = algo.get_results();

        // Eq 6.9 - sum of objective squared must be 1 and  all variables in x_M must be close to 0.5.
        let mut invalid_individuals: usize = 0;
        for ind in &results.individuals {
            let obj_sum: f64 = ind
                .get_objective_values()
                .unwrap()
                .iter()
                .map(|v| v.powi(2))
                .sum();
            assert_approx_eq!(f64, obj_sum, 1.0, epsilon = 0.1);

            let vars: Vec<f64> = ((number_variables - k + 1)..=number_variables)
                .map(|i| {
                    ind.get_variable_value(format!("x{i}").as_str())
                        .unwrap()
                        .as_real()
                        .unwrap()
                })
                .collect();

            let strict_range = 0.480..0.520;
            let values_outside = check_value_in_range(&vars, &strict_range);
            if !values_outside.is_empty() {
                invalid_individuals += 1;
            }
        }

        if invalid_individuals > 10 {
            panic!("Found {invalid_individuals} individuals not meeting the ideal solution bounds",);
        }
    }
}