use std::{convert::Infallible, env, fs, io::Write, rc::Rc};

use eoa_lib::{
    bounded::BoundedOVector, comparison::MinimizingOperator, constraints::{stochastic_ranking_evolution_algorithm, ConstrainedEvalFitness, ConstrainedFitnessFunction, ConstraintFunction, LowerThanConstraintFunction}, crossover::{ArithmeticCrossover, BoundedCrossover, BoundedCrossoverStrategy}, evolution::{EvolutionCandidate, EvolutionResult, EvolutionStats}, fitness::FitnessFunction, initializer::{Initializer, RandomInitializer}, multi_objective_evolution::nsga_2, pairing::AdjacentPairing, perturbation::{BoundedPerturbation, BoundedPerturbationStrategy, MutationPerturbation, RandomDistributionPerturbation}, population::{EvaluatedChromosome, EvaluatedPopulation, Population}, replacement::{BestReplacement, GenerationalReplacement}, selection::{BestSelection, TournamentSelection}
};
use nalgebra::{ArrayStorage, Const, Matrix, SVector};
use rand::RngCore;
use rand_distr::Normal;
use chrono::prelude::*;

pub struct ArbitraryFitness<const SIZE: usize> {
    fun: Box<dyn Fn(SVector<f64, SIZE>) -> f64>
}

impl<const SIZE: usize> ArbitraryFitness<SIZE> {
    pub fn new(fun: Box<dyn Fn(SVector<f64, SIZE>) -> f64>) -> Self {
        Self {
            fun
        }
    }

    pub fn zero() -> Self {
        Self {
            fun: Box::new(|_| 0.0)
        }
    }
}

impl<const SIZE: usize> FitnessFunction for ArbitraryFitness<SIZE> {
    type In = SVector<f64, SIZE>;
    type Out = f64;
    type Err = Infallible;

    fn fit(&self, inp: &Self::In) -> Result<Self::Out, Self::Err> {
        Ok((self.fun)(*inp))
    }
}

/// A constrained optimization problem with clear field definitions
pub struct ConstrainedProblem<const DIM: usize, const CONSTRAINTS: usize> {
    pub name: String,
    pub objective: ArbitraryFitness<DIM>,
    pub constraints: [LowerThanConstraintFunction<SVector<f64, DIM>, f64>; CONSTRAINTS],
    pub bounds: (SVector<f64, DIM>, SVector<f64, DIM>), // (min, max)
    pub optimal_value: f64,
    pub instantiate_fn: Option<Rc<dyn Fn() -> ConstrainedProblem<DIM, CONSTRAINTS>>>
}

impl<const DIM: usize, const CONSTRAINTS: usize> ConstrainedProblem<DIM, CONSTRAINTS> {
    pub fn new(instantiate: Rc<dyn Fn() -> ConstrainedProblem<DIM, CONSTRAINTS>>) -> Self {
        let mut problem = instantiate();
        problem.instantiate_fn = Some(instantiate);
        problem
    }

    pub fn clone(&self) -> Self {
        Self::new(self.instantiate_fn.clone().unwrap())
    }
}

/// Configuration for stochastic ranking method
pub struct StochasticRankingConfig {
    pub population_size: usize,
    pub parents_count: usize,
    pub iterations: usize,
    pub n_param: usize,  // N parameter for stochastic ranking
    pub p_param: f64,    // p parameter for stochastic ranking
    pub mutation_std_dev: f64,
}

/// Configuration for NSGA-II method
pub struct NsgaConfig {
    pub population_size: usize,
    pub parents_count: usize,
    pub iterations: usize,
    pub mutation_std_dev: f64,
}

fn problem_g04() -> ConstrainedProblem<5, 6> {
    ConstrainedProblem::new(Rc::new(||
        ConstrainedProblem {
            name: "g04".to_string(),
            objective: ArbitraryFitness::new(
                Box::new(|vec| {
                    5.3578547 * vec[2].powi(2) + 0.8356891 * vec[0] * vec[4]
                    + 37.293239 * vec[0] - 40792.141
                })
            ),
            constraints: [
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        85.334407 + 0.0056858 * vec[1] * vec[4] + 0.0006262 * vec[0] * vec[3]
                        - 0.0022053 * vec[2] * vec[4] - 92.0
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -85.334407 - 0.0056858 * vec[1] * vec[4] - 0.0006262 * vec[0] * vec[3]
                        + 0.0022053 * vec[2] * vec[4]
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        80.51249 + 0.0071317 * vec[1] * vec[4] + 0.0029955 * vec[0] * vec[1]
                        + 0.0021813 * vec[2].powi(2) - 110.0
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -80.51249 - 0.0071317 * vec[1] * vec[4] - 0.0029955 * vec[0] * vec[1]
                        - 0.0021813 * vec[2].powi(2) + 90.0
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        9.300961 + 0.0047026 * vec[2] * vec[4] + 0.0012547 * vec[0] * vec[2]
                        + 0.0019085 * vec[2] * vec[3] - 25.0
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -9.300961 - 0.0047026 * vec[2] * vec[4] - 0.0012547 * vec[0] * vec[2]
                        - 0.0019085 * vec[2] * vec[3] + 20.0
                    })
                ),
            ],
            bounds: (
                SVector::<f64, 5>::from([78.0, 33.0, 27.0, 27.0, 27.0]),     // min bounds
                SVector::<f64, 5>::from([102.0, 45.0, 45.0, 45.0, 45.0])      // max bounds
            ),
            optimal_value: -30665.53867178333,
            instantiate_fn: None
        }))
}

fn problem_g05() -> ConstrainedProblem<4, 5> {
    ConstrainedProblem::new(Rc::new(||
        ConstrainedProblem {
            name: "g05".to_string(),
            objective: ArbitraryFitness::new(
                Box::new(|vec| {
                    3.0 * vec[0] + 0.000001 * vec[0].powi(3) + 2.0 * vec[1]
                    + (0.000002 / 3.0) * vec[1].powi(3)
                })
            ),
            constraints: [
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -vec[3] + vec[2] - 0.55
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -vec[2] + vec[3] - 0.55
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        1000.0 * (-vec[2] - 0.25).sin() + 1000.0 * (-vec[3] - 0.25) - vec[0] + 894.8
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        1000.0 * (vec[2] - 0.25).sin() + 1000.0 * (vec[2] - vec[3] - 0.25) - vec[1] + 894.8
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        1000.0 * (vec[2] - 0.25).sin() + 1000.0 * (vec[3] - vec[2] - 0.25) + 1294.8
                    })
                ),
            ],
            bounds: (
                SVector::<f64, 4>::from([0.0, 0.0, -0.55, -0.55]),     // min bounds
                SVector::<f64, 4>::from([1200.0, 1200.0, 0.55, 0.55])  // max bounds
            ),
            optimal_value: 5126.4967140071,
            instantiate_fn: None
        }))
}

fn problem_g06() -> ConstrainedProblem<2, 2> {
    ConstrainedProblem::new(Rc::new(||
        ConstrainedProblem {
            name: "g06".to_string(),
            objective: ArbitraryFitness::new(
                Box::new(|vec| (vec[0] - 10.0).powi(3) + (vec[1] - 20.0).powi(3))
            ),
            constraints: [
                LowerThanConstraintFunction::new(
                    Box::new(|vec| -(vec[0] - 5.0).powi(2) - (vec[1] - 5.0).powi(2) + 100.0)
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| (vec[0] - 6.0).powi(2) + (vec[1] - 5.0).powi(2) - 82.81)
                ),
            ],
            bounds: (
                SVector::<f64, 2>::new(0.0, 0.0),     // min bounds
                SVector::<f64, 2>::new(50.0, 50.0)  // max bounds
            ),
            optimal_value: -6961.8137558015,
            instantiate_fn: None
        }))
}

fn problem_g08() -> ConstrainedProblem<2, 2> {
    ConstrainedProblem::new(Rc::new(|| ConstrainedProblem {
        name: "g08".to_string(),
        objective: ArbitraryFitness::new(
            Box::new(|vec| {
                let num = (2.0 * std::f64::consts::PI * vec[0]).sin().powi(3)
                    * (2.0 * std::f64::consts::PI * vec[1]).sin();
                let den = vec[0].powi(3) * (vec[0] + vec[1]);
                -num / den
            })
        ),
        constraints: [
            LowerThanConstraintFunction::new(
                Box::new(move |vec| {
                    let x1 = vec[0];
                    let x2 = vec[1];
                    x1.powi(2) - x2 + 1.0
                })
            ),
            LowerThanConstraintFunction::new(
                Box::new(move |vec| {
                    let x1 = vec[0];
                    let x2 = vec[1];
                    1.0 - x1 + (x2 - 4.0).powi(2)
                })
            ),
        ],
        bounds: (
            SVector::<f64, 2>::new(0.0, 0.0),   // min bounds
            SVector::<f64, 2>::new(10.0, 10.0)  // max bounds
        ),
        optimal_value: -0.0958250414180359,
        instantiate_fn: None
    }))
}

fn problem_g09() -> ConstrainedProblem<7, 4> {
    ConstrainedProblem::new(Rc::new(||
        ConstrainedProblem {
            name: "g09".to_string(),
            objective: ArbitraryFitness::new(
                Box::new(|vec| {
                    (vec[0] - 10.0).powi(2) + 5.0 * (vec[1] - 12.0).powi(2)
                    + vec[2].powi(4) + 3.0 * (vec[3] - 11.0).powi(2)
                    + 10.0 * vec[4].powi(6) + 7.0 * vec[5].powi(2)
                    + vec[6].powi(4) - 4.0 * vec[5] * vec[6]
                    - 10.0 * vec[5] - 8.0 * vec[6]
                })
            ),
            constraints: [
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -127.0 + 2.0 * vec[0].powi(2) + 3.0 * vec[1].powi(4)
                        + vec[2] + 4.0 * vec[3].powi(2) + 5.0 * vec[4]
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -282.0 + 7.0 * vec[0] + 3.0 * vec[1] + 10.0 * vec[2].powi(2)
                        + vec[3] - vec[4]
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        -196.0 + 23.0 * vec[0] + vec[1].powi(2) + 6.0 * vec[5].powi(2)
                        - 8.0 * vec[6]
                    })
                ),
                LowerThanConstraintFunction::new(
                    Box::new(|vec| {
                        4.0 * vec[0].powi(2) + vec[1].powi(2) - 3.0 * vec[0] * vec[1]
                        + 2.0 * vec[2].powi(2) + 5.0 * vec[5] - 11.0 * vec[6]
                    })
                ),
            ],
            bounds: (
                SVector::<f64, 7>::from([-10.0; 7]),     // min bounds (all -10)
                SVector::<f64, 7>::from([10.0; 7])       // max bounds (all 10)
            ),
            optimal_value: 680.6300573745,
            instantiate_fn: None
        }))
}

pub fn problem_g11(eps: f64) -> ConstrainedProblem<2, 1> {
    let problem = ConstrainedProblem::new(Rc::new(move || ConstrainedProblem {
        name: "g11".to_string(),
        objective: ArbitraryFitness::new(
            Box::new(|vec| {
                // Minimize f(x) = x1^2 + (x2 - 1)^2
                vec[0].powi(2) + (vec[1] - 1.0).powi(2)
            })
        ),
        constraints: [
            // Equality h(x) = x2 - x1^2 = 0
            // |h| - eps >= 0.0
            LowerThanConstraintFunction::new(
                Box::new(move |vec| {
                    let h = vec[1] - vec[0].powi(2);
                    h.abs() - eps
                })
            ),
        ],
        bounds: (
            SVector::<f64, 2>::new(-50.0, -50.0),  // min bounds
            SVector::<f64, 2>::new(50.0, 50.0)     // max bounds
        ),
        optimal_value: 0.7499, // Best known optimum
        instantiate_fn: None
    }));

    problem
}

fn problem_g21() -> ConstrainedProblem<7, 6> {
    ConstrainedProblem::new(Rc::new(|| ConstrainedProblem {
        name: "g21".to_string(),
        objective: ArbitraryFitness::new(Box::new(|vec| vec[0])),
        constraints: [
            LowerThanConstraintFunction::new(Box::new(|vec| {
                -vec[0] + 35.0 * vec[1].powf(0.6) + 35.0 * vec[2].powf(0.6)
            })),
            LowerThanConstraintFunction::new(Box::new(|vec| {
                -300.0 * vec[2] + 7500.0 * vec[4] - 7500.0 * vec[5]
                - 25.0 * vec[3] * vec[4] + 25.0 * vec[3] * vec[5] + vec[2] * vec[3]
            })),
            LowerThanConstraintFunction::new(Box::new(|vec| {
                100.0 * vec[1] + 155.365 * vec[3] + 2500.0 * vec[6]
                - vec[1] * vec[3] - 25.0 * vec[3] * vec[6] - 15536.5
            })),
            LowerThanConstraintFunction::new(Box::new(|vec| {
                -vec[4] + (-vec[3] + 900.0).ln()
            })),
            LowerThanConstraintFunction::new(Box::new(|vec| {
                 -vec[5] + (vec[3] + 300.0).ln()
            })),
            LowerThanConstraintFunction::new(Box::new(|vec| {
                 -vec[6] + (- 2.0 * vec[3] + 700.0).ln()
            })),
        ],
        bounds: (
            // x1..x7 Min
            SVector::<f64, 7>::from([0.0, 0.0, 0.0, 100.0, 6.3, 5.9, 4.5]),
            // x1..x7 Max
            SVector::<f64, 7>::from([1000.0, 40.0, 40.0, 300.0, 6.7, 6.4, 6.25])
        ),
        optimal_value: 193.724510070035,
        instantiate_fn: None
    }))
}

pub fn problem_g24() -> ConstrainedProblem<2, 2> {
    ConstrainedProblem::new(Rc::new(|| ConstrainedProblem {
        name: "g24".to_string(),
        objective: ArbitraryFitness::new(
            Box::new(|vec| {
                // Minimize f(x) = -x1 - x2
                -vec[0] - vec[1]
            })
        ),
        constraints: [
            // g1(x) = -2x1^4 + 8x1^3 - 8x1^2 + x2 - 2 <= 0
            LowerThanConstraintFunction::new(
                Box::new(|vec| {
                    -2.0 * vec[0].powi(4) + 8.0 * vec[0].powi(3) - 8.0 * vec[0].powi(2) + vec[1] - 2.0
                })
            ),
            // g2(x) = -4x1^4 + 32x1^3 - 88x1^2 + 96x1 + x2 - 36 <= 0
            LowerThanConstraintFunction::new(
                Box::new(|vec| {
                    -4.0 * vec[0].powi(4) + 32.0 * vec[0].powi(3) - 88.0 * vec[0].powi(2) + 96.0 * vec[0] + vec[1] - 36.0
                })
            ),
        ],
        bounds: (
            SVector::<f64, 2>::new(0.0, 0.0),  // min bounds
            SVector::<f64, 2>::new(3.0, 4.0)   // max bounds
        ),
        optimal_value: -5.50801327159536, // Best known optimum
        instantiate_fn: None
    }))
}

/// Solve a constrained optimization problem using stochastic ranking
///
/// Returns the evolution result with feasible fractions for each iteration
pub fn solve_with_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>(
    mut problem: ConstrainedProblem<DIM, CONSTRAINTS>,
    population_size: usize,
    parents_count: usize,
    iterations: usize,
    stochastic_params: (usize, f64), // (N, p) for stochastic ranking
    mutation_std_dev: f64,
    rng: &mut dyn RngCore,
) -> Result<(EvolutionResult<SVector<f64, DIM>, f64>, Vec<f64>, Vec<f64>), Box<dyn std::error::Error>> {
    // Create initial population
    let initializer = RandomInitializer::new(
        Box::new(BoundedOVector::new(problem.bounds.0, problem.bounds.1))
    );
    let initial_population = Population::from_vec(
        initializer.initialize(nalgebra::Const::<DIM>, population_size, rng)
    );

    // Setup components as specified
    // let mut selection = TournamentSelection::new(5, 0.9);
    let mut selection = TournamentSelection::new(5, 0.95);
    let mut replacement = GenerationalReplacement;
    let mut pairing = AdjacentPairing::new();
    let mut crossover = ArithmeticCrossover::new();
    // let mut crossover = BoundedCrossover::<nalgebra::Const<2>, 2, _>::new(
    //     ArithmeticCrossover::new(),
    //     problem.bounds.0,
    //     problem.bounds.1,
    //     BoundedCrossoverStrategy::Retry(5)
    // );

    // Setup bounded random distribution perturbation with Normal distribution
    let normal_perturbation = RandomDistributionPerturbation::<DIM, Normal<f64>>::normal(mutation_std_dev)?;
    let mut perturbation = normal_perturbation;
    // let perturbation = BoundedPerturbation::new(
    //     normal_perturbation,
    //     problem.bounds.0,
    //     problem.bounds.1,
    //     BoundedPerturbationStrategy::Retry(5)
    // );
    let mut mutation = MutationPerturbation::new(Box::new(perturbation), 0.1);

    // The weight is so large mainly because of the g11 that has very small values.
    // Somehow the higher weights do seem to help, even though I am unsure why exactly.
    let constraint_weights = [1.0; CONSTRAINTS];

    let (N, p) = stochastic_params;
    let better_than = MinimizingOperator::new();

    // Convert constraint array references
    let constraint_refs = problem.constraints.iter().collect::<Vec<_>>().try_into()
        .map_err(|_| "Failed to convert constraint references")?;

    let result = stochastic_ranking_evolution_algorithm(
        initial_population,
        parents_count,
        N,
        p,
        &mut problem.objective,
        constraint_refs,
        constraint_weights,
        &mut pairing,
        &mut selection,
        &mut crossover,
        &mut mutation,
        &mut replacement,
        &better_than,
        iterations,
        rng)?;

    // Extract feasible fractions from the result
    let (evolution_result, feasible_fractions) = result;
    
    // For now, create placeholder constraint violations data
    // TODO: This needs library-level changes to properly track constraint violations
    let avg_constraint_violations = vec![0.0; feasible_fractions.len()];
    
    Ok((evolution_result, feasible_fractions, avg_constraint_violations))
}

/// Helper function to check if a chromosome is feasible
fn check_feasibility<const DIM: usize, const CONSTRAINTS: usize>(
    chromosome: &SVector<f64, DIM>,
    constraints: &[LowerThanConstraintFunction<SVector<f64, DIM>, f64>; CONSTRAINTS],
) -> bool {
    constraints.iter().all(|constraint| {
        constraint.is_feasible(chromosome).unwrap_or(false)
    })
}

/// Individual constraint fitness function - wraps a single constraint to act as a fitness function
pub struct SingleConstraintFitness<const DIM: usize, TConstraint: ConstraintFunction<Chromosome = SVector<f64, DIM>, Out = f64>> {
    constraint: TConstraint,
    capped: bool
}

impl<const DIM: usize, TConstraint: ConstraintFunction<Chromosome = SVector<f64, DIM>, Out = f64>> SingleConstraintFitness<DIM, TConstraint> {
    pub fn new(constraint: TConstraint, capped: bool) -> Self {
        Self { constraint, capped }
    }
}

impl<const DIM: usize, TConstraint: ConstraintFunction<Chromosome = SVector<f64, DIM>, Out = f64>> FitnessFunction for SingleConstraintFitness<DIM, TConstraint> {
    type In = SVector<f64, DIM>;
    type Out = f64;
    type Err = Infallible;

    fn fit(&self, inp: &Self::In) -> Result<Self::Out, Self::Err> {
        Ok(if self.constraint.is_feasible(inp).unwrap() && self.capped {
            0.0
        } else {
            self.constraint.evaluate(inp).unwrap()
        })
    }
}

/// Solve a constrained optimization problem using NSGA-II
pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
    problem: ConstrainedProblem<DIM, CONSTRAINTS>,
    population_size: usize,
    parents_count: usize,
    iterations: usize,
    mutation_std_dev: f64,
    rng: &mut dyn RngCore,
) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; 2]>, Vec<f64>, Vec<f64>), Box<dyn std::error::Error>> {
    // Create initial population
    let initializer = RandomInitializer::new(
        Box::new(BoundedOVector::new(problem.bounds.0, problem.bounds.1))
    );
    let initial_population = Population::from_vec(
        initializer.initialize(nalgebra::Const::<DIM>, population_size, rng)
    );

    // Setup components
    let mut pairing = AdjacentPairing::new();
    let mut crossover = ArithmeticCrossover::new();

    // Setup bounded random distribution perturbation with Normal distribution
    let normal_perturbation = RandomDistributionPerturbation::<DIM, Normal<f64>>::normal(mutation_std_dev)?;
    let mut mutation = MutationPerturbation::new(Box::new(normal_perturbation), 0.1);

    let better_than = MinimizingOperator::new();

    // Create objectives: both using ConstrainedFitnessFunction with different fitness components
    let constraint_weights = [1.0e9; CONSTRAINTS];
    let constraint_refs = problem.constraints.iter().collect::<Vec<_>>().try_into()
        .map_err(|_| "Failed to convert constraint references")?;

    let zero_fitness = ArbitraryFitness::zero();

    // Second objective: constraint violation only (zero fitness + constraints)
    let constrained_fitness_obj2 = ConstrainedFitnessFunction {
        fitness: &zero_fitness,
        constraints: constraint_refs,
        constraint_weights,
        capped: true
    };

    let objectives: [Box<dyn FitnessFunction<In = SVector<f64, DIM>, Out = f64, Err = Infallible>>; 2] = [
        Box::new(problem.objective),
        Box::new(constrained_fitness_obj2),
    ];

    let mut feasible_fractions = Vec::with_capacity(iterations);
    let mut avg_constraint_violations = Vec::with_capacity(iterations);

    let result = nsga_2(
        initial_population,
        parents_count,
        objectives,
        &mut pairing,
        &mut crossover,
        &mut mutation,
        &better_than,
        iterations,
        rng,
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
            // Calculate feasible fraction based on second objective being close to zero
            let feasible_count = population.population
                .iter()
                .filter(|individual| {
                    individual.evaluation.evaluations[1] <= 0.0
                })
                .count();

            // Calculate average constraint violation (second objective is constraint violation)
            let avg_constraint_violation = population.population
                .iter()
                .map(|individual| {
                    individual.evaluation.evaluations[1].max(0.0) // Only positive values are violations
                })
                .sum::<f64>() / population.population.len() as f64;

            let feasible_fraction = feasible_count as f64 / population.population.len() as f64;
            feasible_fractions.push(feasible_fraction);
            avg_constraint_violations.push(avg_constraint_violation);
        },
        |_, evaluation, best_candidate| {
            // Do not save infeasible solutions!
            if evaluation.evaluations[1] > 0.0 {
                return false;
            }

            if best_candidate.is_none() {
                return true;
            }

            evaluation.evaluations[0] < best_candidate.as_ref().unwrap().evaluated_chromosome.evaluation.evaluations[0]
        }
    )?;

    Ok((result, feasible_fractions, avg_constraint_violations))
}

/// Solve a constrained optimization problem using NSGA-II with individual constraint objectives
/// For simplicity, this function only works with 2-constraint problems
pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const CONSTRS_PLUS_ONE: usize>(
    problem: ConstrainedProblem<DIM, CONSTRAINTS>,
    population_size: usize,
    parents_count: usize,
    iterations: usize,
    mutation_std_dev: f64,
    rng: &mut dyn RngCore,
    capped: bool
) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>, Vec<f64>, Vec<f64>), Box<dyn std::error::Error>> {
    // Unfortunately Rustc doesn't support addition in generics...
    assert_eq!(CONSTRAINTS + 1, CONSTRS_PLUS_ONE);

    // Clone the problem to get bounds info first
    let bounds = (problem.bounds.0, problem.bounds.1);

    // Create initial population
    let initializer = RandomInitializer::new(
        Box::new(BoundedOVector::new(bounds.0, bounds.1))
    );
    let initial_population = Population::from_vec(
        initializer.initialize(nalgebra::Const::<DIM>, population_size, rng)
    );

    // Setup components
    let mut pairing = AdjacentPairing::new();
    let mut crossover = ArithmeticCrossover::new();
    let normal_perturbation = RandomDistributionPerturbation::<DIM, Normal<f64>>::normal(mutation_std_dev)?;
    let mut mutation = MutationPerturbation::new(Box::new(normal_perturbation), 0.1);
    let better_than = MinimizingOperator::new();

    // Create objectives: fitness + individual constraints using cloned problem
    let mut objective = Some(problem.objective);
    let mut constraints = problem.constraints.into_iter();

    let objectives: [Box<dyn FitnessFunction<In = SVector<f64, DIM>, Out = f64, Err = Infallible>>; CONSTRS_PLUS_ONE] = std::array::from_fn(move |i| {
        let val: Box<dyn FitnessFunction<In = SVector<f64, DIM>, Out = f64, Err = Infallible>> = if i == 0 {
            let obj = objective.take().expect("Taken already!");
            Box::new(obj)
        } else {
            Box::new(SingleConstraintFitness::new(constraints.next().unwrap(), capped))
        };

        val
    });

    let mut feasible_fractions = Vec::with_capacity(iterations);
    let mut avg_constraint_violations = Vec::with_capacity(iterations);

    let result = nsga_2(
        initial_population,
        parents_count,
        objectives,
        &mut pairing,
        &mut crossover,
        &mut mutation,
        &better_than,
        iterations,
        rng,
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
            // Calculate feasible fraction - all constraints (objectives 1 and 2) must be <= 0
            let feasible_count: f64 =
                population.population.iter().filter(
                    |individual| {
                        individual.evaluation.evaluations
                            .iter()
                            .skip(1)
                            .all(|&eval| eval <= 0.0)
                    }
                ).count() as f64;

            // Calculate average constraint violation (sum all constraint violations from objectives 1+)
            let total_violation: f64 = population.population
                .iter()
                .map(|individual| {
                    individual.evaluation.evaluations
                        .iter()
                        .skip(1) // Skip fitness objective, only look at constraints
                        .map(|&eval| eval.max(0.0)) // Only positive values are violations
                        .sum::<f64>()
                })
                .sum();
            let avg_constraint_violation = total_violation / population.population.len() as f64;

            let feasible_fraction = feasible_count as f64 / population.population.len() as f64;
            feasible_fractions.push(feasible_fraction);
            avg_constraint_violations.push(avg_constraint_violation);
        },
        |_, evaluation, best_candidate| {
            // Only save feasible solutions (all constraints satisfied)
            // Skip the first objective (fitness), check constraints (objectives 1+)
            if evaluation.evaluations.iter().skip(1).any(|&eval| eval > 0.0) {
                return false;
            }

            if best_candidate.is_none() {
                return true;
            }

            // Compare based on fitness (first objective)
            evaluation.evaluations[0] < best_candidate.as_ref().unwrap().evaluated_chromosome.evaluation.evaluations[0]
        }
    )?;

    Ok((result, feasible_fractions, avg_constraint_violations))
}

const ITERATIONS: usize = 1000;
const POPULATION: usize = 500;
const PARENTS_COUNT: usize = 500;
const G11_EPS: f64 = 0.00015;

fn handle_g06_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g06();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 1.0,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g08_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g08();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 0.5,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g11_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g11(G11_EPS);
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: POPULATION * 2,

        p_param: 0.45,
        mutation_std_dev: 0.01 / 50.0,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g04_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g04();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 1.0,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g05_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g05();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 10.0,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g09_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g09();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 1.0,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g21_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g21();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 10.0,
    };
    run_stochastic_ranking(problem, config)
}

fn handle_g24_srank() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g24();
    let config = StochasticRankingConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        n_param: 2 * POPULATION,
        p_param: 0.45,
        mutation_std_dev: 0.1,
    };
    run_stochastic_ranking(problem, config)
}

/// Generic function to save evolution results to CSV files with date-based naming
fn save_evolution_results<const DIM: usize, const CONSTRAINTS: usize, TEval: std::fmt::Debug>(
    method: &str,
    problem: &ConstrainedProblem<DIM, CONSTRAINTS>,
    evolution_result: &EvolutionResult<SVector<f64, DIM>, TEval>,
    feasible_fractions: &[f64],
    avg_constraint_violations: &[f64],
) -> Result<(), Box<dyn std::error::Error>> {
    // Get current date and time for unique file naming
    let now = Local::now();
    let timestamp = now.format("%Y%m%d_%H%M%S").to_string();

    // Create output directories
    let output_dir = format!("solutions/{}/{}", method, problem.name);
    let feasible_dir = format!("{}/feasible_fraction", output_dir);
    let constraint_dir = format!("{}/constraint_violation", output_dir);
    fs::create_dir_all(&output_dir)?;
    fs::create_dir_all(&feasible_dir)?;
    fs::create_dir_all(&constraint_dir)?;

    // Write best candidates CSV with timestamp
    let best_candidates_path = format!("{}/best_candidates_{}.csv", output_dir, timestamp);
    let mut best_file = fs::File::create(&best_candidates_path)?;
    writeln!(best_file, "iteration,evaluation,fitness")?;

    // Write evolution stats (best candidates through iterations)
    for candidate in &evolution_result.stats.best_candidates {
        writeln!(best_file, "{},{:?}", candidate.evaluation, candidate.evaluated_chromosome.evaluation)?;
    }

    // Write final best candidate and total evaluations
    if let Some(ref best_candidate) = evolution_result.best_candidate {
        writeln!(best_file, "{},{:?}", evolution_result.evaluations, best_candidate.evaluation)?;
    }

    // Write feasible fractions CSV with timestamp
    let feasible_path = format!("{}/feasible_fractions_{}.csv", feasible_dir, timestamp);
    let mut feasible_file = fs::File::create(&feasible_path)?;
    writeln!(feasible_file, "iteration,feasible_fraction")?;

    for (i, fraction) in feasible_fractions.iter().enumerate() {
        writeln!(feasible_file, "{},{}", i, fraction)?;
    }

    // Write constraint violations CSV with timestamp
    let constraint_path = format!("{}/constraint_violations_{}.csv", constraint_dir, timestamp);
    let mut constraint_file = fs::File::create(&constraint_path)?;
    writeln!(constraint_file, "iteration,avg_constraint_violation")?;

    for (i, violation) in avg_constraint_violations.iter().enumerate() {
        writeln!(constraint_file, "{},{}", i, violation)?;
    }

    println!("Results saved to:");
    println!("  Best candidates: {}", best_candidates_path);
    println!("  Feasible fractions: {}", feasible_path);
    println!("  Constraint violations: {}", constraint_path);

    Ok(())
}

fn run_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>(
    problem: ConstrainedProblem<DIM, CONSTRAINTS>,
    config: StochasticRankingConfig,
) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();

    let result = solve_with_stochastic_ranking(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        (config.n_param, config.p_param),
        config.mutation_std_dev,
        &mut rng,
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;

    // Save results to CSV files
    save_evolution_results("srank", &problem, &evolution_result, &feasible_fractions, &avg_constraint_violations)?;

    if let Some(best_candidate) = evolution_result.best_candidate {
        println!("Best solution found:");
        println!("  Chromosome: {:?}", best_candidate.chromosome);
        println!(
            "  Fitness: {} ({} %)",
            best_candidate.evaluation,
            ((best_candidate.evaluation - problem.optimal_value) / problem.optimal_value).abs() * 100.0);
        println!("  Iterations: {}", evolution_result.iterations);
        println!("  Evaluations: {}", evolution_result.evaluations);
        println!("  Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);

        // Sanity check: verify the best candidate is feasible
        let best_chromosome = &best_candidate.chromosome;

        println!("\nFeasibility check for best solution:");
        for (i, constraint) in problem.constraints.iter().enumerate() {
            match constraint.evaluate(best_chromosome) {
                Ok(value) => {
                    let is_feasible = value <= 0.0;
                    println!("  Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
                }
                Err(e) => {
                    println!("  Constraint {}: Error evaluating - {}", i+1, e);
                }
            }
        }
    } else {
        println!("Could not find any feasible solution!")
    }

    Ok(())
}

fn run_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
    problem: ConstrainedProblem<DIM, CONSTRAINTS>,
    config: NsgaConfig,
) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();

    let result = solve_with_nsga_ii(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;

    // Save results to CSV files
    save_evolution_results("nsga", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?;

    if let Some(best_candidate) = evolution_result.best_candidate {
        println!("Best solution found:");
        println!("  Chromosome: {:?}", best_candidate.chromosome);
        println!("  Objectives: {:?} ({} %)", best_candidate.evaluation, ((best_candidate.evaluation[0] - problem.optimal_value) / problem.optimal_value).abs() * 100.0);
        println!("  Iterations: {}", evolution_result.iterations);
        println!("  Evaluations: {}", evolution_result.evaluations);
        println!("  Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);

        // Sanity check: verify the best candidate is feasible
        let best_chromosome = &best_candidate.chromosome;

        println!("\nFeasibility check for best solution:");
        for (i, constraint) in problem.constraints.iter().enumerate() {
            match constraint.evaluate(best_chromosome) {
                Ok(value) => {
                    let is_feasible = value <= 0.0;
                    println!("  Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
                }
                Err(e) => {
                    println!("  Constraint {}: Error evaluating - {}", i+1, e);
                }
            }
        }
    } else {
        println!("Could not find any feasible solution!")
    }

    Ok(())
}

fn run_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const CONSTRS_PLUS_ONE: usize>(
    problem: ConstrainedProblem<DIM, CONSTRAINTS>,
    evolution_result: EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>,
    feasible_fractions: Vec<f64>,
    avg_constraint_violations: Vec<f64>,
    capped: bool
) -> Result<(), Box<dyn std::error::Error>> {
    // Unfortunately Rustc doesn't support addition in generics...
    assert_eq!(CONSTRAINTS + 1, CONSTRS_PLUS_ONE);

    // Save results to CSV files
    if capped {
        save_evolution_results("nsga_multi", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?;
    } else {
        save_evolution_results("nsga_multi_noncapped", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?;
    }

    if let Some(best_candidate) = evolution_result.best_candidate {
        println!("Best solution found:");
        println!("  Chromosome: {:?}", best_candidate.chromosome);
        println!("  Objectives: {:?} ({} %)", best_candidate.evaluation, ((best_candidate.evaluation[0] - problem.optimal_value) / problem.optimal_value).abs() * 100.0);
        println!("  Iterations: {}", evolution_result.iterations);
        println!("  Evaluations: {}", evolution_result.evaluations);
        println!("  Final feasible fraction: {:.2}%", feasible_fractions.last().unwrap_or(&0.0) * 100.0);

        // Sanity check: verify the best candidate is feasible
        let best_chromosome = &best_candidate.chromosome;

        println!("\nFeasibility check for best solution:");
        for (i, constraint) in problem.constraints.iter().enumerate() {
            match constraint.evaluate(best_chromosome) {
                Ok(value) => {
                    let is_feasible = value <= 0.0;
                    println!("  Constraint {}: {} ({})", i+1, value, if is_feasible { "FEASIBLE" } else { "INFEASIBLE" });
                }
                Err(e) => {
                    println!("  Constraint {}: Error evaluating - {}", i+1, e);
                }
            }
        }
    } else {
        println!("Could not find any feasible solution!")
    }

    Ok(())
}

// NSGA-II handler functions
fn handle_g06_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g06();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 1.0,
    };
    run_nsga_ii(problem, config)
}

fn handle_g08_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g08();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.5,
    };
    run_nsga_ii(problem, config)
}

fn handle_g11_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g11(G11_EPS);
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.01 / 50.0,
    };
    run_nsga_ii(problem, config)
}

fn handle_g04_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g04();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 1.0,
    };
    run_nsga_ii(problem, config)
}

fn handle_g05_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g05();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 10.0,
    };
    run_nsga_ii(problem, config)
}

fn handle_g09_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g09();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 1.0,
    };
    run_nsga_ii(problem, config)
}

fn handle_g21_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g21();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 10.0,
    };
    run_nsga_ii(problem, config)
}

fn handle_g24_nsga() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g24();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.1,
    };
    run_nsga_ii(problem, config)
}

// NSGA-Multi handler functions for individual problems
fn handle_nsga_multi_g06(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g06();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.5,
    };

    let result = solve_with_nsga_multi::<2, 2, 3>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g08(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g08();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.5,
    };

    let result = solve_with_nsga_multi::<2, 2, 3>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g11(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g11(G11_EPS);
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.01 / 50.0,
    };

    let result = solve_with_nsga_multi::<2, 1, 2>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g04(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g04();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 1.0,
    };

    let result = solve_with_nsga_multi::<5, 6, 7>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g05(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g05();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 10.0,
    };

    let result = solve_with_nsga_multi::<4, 5, 6>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g09(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g09();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 1.0,
    };

    let result = solve_with_nsga_multi::<7, 4, 5>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g21(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g21();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 10.0,
    };

    let result = solve_with_nsga_multi::<7, 6, 7>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn handle_nsga_multi_g24(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
    let mut rng = rand::rng();
    let problem = problem_g24();
    let config = NsgaConfig {
        population_size: POPULATION,
        parents_count: PARENTS_COUNT,
        iterations: ITERATIONS,
        mutation_std_dev: 0.1,
    };

    let result = solve_with_nsga_multi::<2, 2, 3>(
        problem.clone(),
        config.population_size,
        config.parents_count,
        config.iterations,
        config.mutation_std_dev,
        &mut rng,
        capped
    )?;

    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
}

fn main() {
    let args: Vec<String> = env::args().collect();

    if args.len() != 3 {
        eprintln!("Usage: {} <method> <problem>", args[0]);
        eprintln!("Methods: srank, nsga, nsga_multi, nsga_multi_noncapped");
        eprintln!("Problems: g04, g05, g06, g08, g09, g11, g21, g24");
        std::process::exit(1);
    }

    let method = &args[1];
    let problem = &args[2];

    let result = match (method.as_str(), problem.as_str()) {
        ("srank", "g04") => handle_g04_srank(),
        ("srank", "g05") => handle_g05_srank(),
        ("srank", "g06") => handle_g06_srank(),
        ("srank", "g08") => handle_g08_srank(),
        ("srank", "g09") => handle_g09_srank(),
        ("srank", "g11") => handle_g11_srank(),
        ("srank", "g21") => handle_g21_srank(),
        ("srank", "g24") => handle_g24_srank(),
        ("nsga", "g04") => handle_g04_nsga(),
        ("nsga", "g05") => handle_g05_nsga(),
        ("nsga", "g06") => handle_g06_nsga(),
        ("nsga", "g08") => handle_g08_nsga(),
        ("nsga", "g09") => handle_g09_nsga(),
        ("nsga", "g11") => handle_g11_nsga(),
        ("nsga", "g21") => handle_g21_nsga(),
        ("nsga", "g24") => handle_g24_nsga(),
        ("nsga_multi", "g04") => handle_nsga_multi_g04(true),
        ("nsga_multi", "g05") => handle_nsga_multi_g05(true),
        ("nsga_multi", "g06") => handle_nsga_multi_g06(true),
        ("nsga_multi", "g08") => handle_nsga_multi_g08(true),
        ("nsga_multi", "g09") => handle_nsga_multi_g09(true),
        ("nsga_multi", "g11") => handle_nsga_multi_g11(true),
        ("nsga_multi", "g21") => handle_nsga_multi_g21(true),
        ("nsga_multi", "g24") => handle_nsga_multi_g24(true),
        ("nsga_multi_noncapped", "g04") => handle_nsga_multi_g04(false),
        ("nsga_multi_noncapped", "g05") => handle_nsga_multi_g05(false),
        ("nsga_multi_noncapped", "g06") => handle_nsga_multi_g06(false),
        ("nsga_multi_noncapped", "g08") => handle_nsga_multi_g08(false),
        ("nsga_multi_noncapped", "g09") => handle_nsga_multi_g09(false),
        ("nsga_multi_noncapped", "g11") => handle_nsga_multi_g11(false),
        ("nsga_multi_noncapped", "g21") => handle_nsga_multi_g21(false),
        ("nsga_multi_noncapped", "g24") => handle_nsga_multi_g24(false),
        (_, _) => {
            eprintln!("Invalid method '{}' or problem '{}'", method, problem);
            eprintln!("Methods: srank, nsga, nsga_multi, nsga_multi_noncapped");
            eprintln!("Problems: g04, g05, g06, g08, g09, g11, g21, g24");
            std::process::exit(1);
        }
    };

    if let Err(e) = result {
        eprintln!("Error: {}", e);
        std::process::exit(1);
    }
}