~ruther/ctu-fee-eoa

db937f4c8c792edfa6e98b20bc44f5036bc64699 — Rutherther 6 days ago f1a6bbd
feat: implement improved nsga for constraints with archive of feasible solutions
2 files changed, 344 insertions(+), 38 deletions(-)

M codes/constr_hw02/src/main.rs
M codes/eoa_lib/src/multi_objective_evolution.rs
M codes/constr_hw02/src/main.rs => codes/constr_hw02/src/main.rs +336 -27
@@ 1,4 1,4 @@
use std::{convert::Infallible, env, fs, io::Write, rc::Rc};
use std::{convert::Infallible, env, fs, io::Write, rc::Rc, cell::RefCell};

use eoa_lib::{
    bounded::BoundedOVector, comparison::MinimizingOperator, constraints::{stochastic_ranking_evolution_algorithm, ConstrainedEvalFitness, ConstrainedFitnessFunction, ConstraintFunction, LowerThanConstraintFunction}, crossover::{ArithmeticCrossover, BoundedCrossover, BoundedCrossoverStrategy, Crossover}, evolution::{EvolutionCandidate, EvolutionResult, EvolutionStats}, fitness::FitnessFunction, initializer::{Initializer, RandomInitializer}, multi_objective_evolution::{constrained_nsga_2, nsga_2, NSGAEvaluation}, pairing::{AdjacentPairing, ParentPairing}, perturbation::{BoundedPerturbation, BoundedPerturbationStrategy, MutationPerturbation, RandomDistributionPerturbation}, population::{EvaluatedChromosome, EvaluatedPopulation, Population}, replacement::{BestReplacement, GenerationalReplacement}, selection::{BestSelection, TournamentSelection}


@@ 117,34 117,37 @@ impl<
            let a_feasible = a.evaluation.evaluations.iter().skip(1).all(|&constr| constr <= 0.0);
            let b_feasible = b.evaluation.evaluations.iter().skip(1).all(|&constr| constr <= 0.0);

            match (a_feasible, b_feasible) {
                (false, true) => {
                    if rng.random_bool(self.p_single_replaced) {
                        pair[0] = rng.random_range(parents_count..full_population);
            // Only proceed with replacement if we have archived feasible solutions
            if full_population > parents_count {
                match (a_feasible, b_feasible) {
                    (false, true) => {
                        if rng.random_bool(self.p_single_replaced) {
                            pair[0] = rng.random_range(parents_count..full_population);
                        }
                    },
                    (true, false) => {
                        if rng.random_bool(self.p_single_replaced) {
                            pair[1] = rng.random_range(parents_count..full_population);
                        }
                    },
                    (false, false) => {
                        if rng.random_bool(self.p_double_first_replaced) {
                            pair[0] = rng.random_range(parents_count..full_population);

                            if rng.random_bool(self.p_double_second_replaced) {
                              pair[1] = rng.random_range(parents_count..full_population);
                            }
                        }
                    }
                },
                (true, false) => {
                    if rng.random_bool(self.p_single_replaced) {
                        pair[1] = rng.random_range(parents_count..full_population);
                    (true, true) => {
                        // Do nothing.
                    }
                },
                (false, false) => {
                    if rng.random_bool(self.p_double_first_replaced) {
                        pair[0] = rng.random_range(parents_count..full_population);

                        if rng.random_bool(self.p_double_second_replaced) {
                          pair[1] = rng.random_range(parents_count..full_population);
                        }
                    }
                }
                (true, true) => {
                    // Do nothing.
                }

            }
        }

        self.crossover(&joined_population, new_pairs.into_iter(), rng)
        self.crossover.crossover(&joined_population, new_pairs.into_iter(), rng)
    }
}



@@ 677,7 680,7 @@ pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
        &better_than,
        iterations,
        rng,
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>, _crossover: &mut _| {
            // Calculate feasible fraction based on second objective being close to zero
            let feasible_count = population.population
                .iter()


@@ 775,7 778,7 @@ pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const C
        &better_than,
        iterations,
        rng,
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>, _crossover: &mut _| {
            // Calculate feasible fraction - all constraints (objectives 1 and 2) must be <= 0
            let feasible_count: f64 =
                population.population.iter().filter(


@@ 893,7 896,7 @@ pub fn solve_with_nsga_constr<const DIM: usize, const CONSTRAINTS: usize>(
        &better_than,
        iterations,
        rng,
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>| {
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>, _crossover: &mut _| {
            // Calculate feasible fraction and average constraint violation
            let feasible_count = population.population
                .iter()


@@ 941,8 944,149 @@ pub fn solve_with_nsga_constr<const DIM: usize, const CONSTRAINTS: usize>(
    Ok((result, feasible_fractions, avg_constraint_violations))
}

pub fn solve_with_nsga_improved<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();

    // Create the wrapped crossover with arithmetic crossover inside
    let mut wrapped_crossover = FeasibleCrossoverWrapper {
        p_single_replaced: P_SINGLE_REPLACED,
        p_double_first_replaced: P_DOUBLE_FIRST_REPLACED,
        p_double_second_replaced: P_DOUBLE_SECOND_REPLACED,
        archived_count: ARCHIVE_SIZE,
        archived_population: Vec::new(),
        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();

    // Clone the problem to access its fields
    let cloned_problem = problem.clone();

    // Create objectives: following the nsga algorithm pattern
    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(cloned_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::<2, SVector<f64, DIM>, f64, Infallible, 2, _, _, _>(
        initial_population,
        parents_count,
        objectives,
        &mut pairing,
        &mut wrapped_crossover,
        &mut mutation,
        &better_than,
        iterations,
        rng,
        |_iteration: usize, _stats: &EvolutionStats<SVector<f64, DIM>, _>, population: &EvaluatedPopulation<SVector<f64, DIM>, _>, crossover: &mut _| {
            // Update archive with feasible solutions
            let feasible_individuals: Vec<_> = population.population.iter()
                .filter(|individual| {
                    problem.constraints.iter().all(|constraint| {
                        constraint.is_feasible(&individual.chromosome).unwrap_or(false)
                    })
                })
                .cloned()
                .collect();

            // Update the crossover's archive now that we have access to it
            crossover.update_archive(population);

            // Calculate feasible fraction and average constraint violation
            let feasible_count = population.population
                .iter()
                .filter(|individual| {
                    // Check if the individual satisfies all constraints using the actual constraint functions
                    problem.constraints.iter().all(|constraint| {
                        constraint.is_feasible(&individual.chromosome).unwrap_or(false)
                    })
                })
                .count();

            // Calculate average constraint violation
            let total_violation: f64 = population.population
                .iter()
                .map(|individual| {
                    problem.constraints
                        .iter()
                        .map(|constraint| constraint.evaluate(&individual.chromosome).unwrap_or(0.0).max(0.0))
                        .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);
        },
        |chromosome, evaluation, best_candidate| {
            // Only save feasible solutions using actual constraint functions
            if !problem.constraints.iter().all(|constraint| {
                constraint.is_feasible(chromosome).unwrap_or(false)
            }) {
                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]
        }
    )?;

    // Return both objectives (original objective and constraint sum)
    let objective_result = result;

    Ok((objective_result, feasible_fractions, avg_constraint_violations))
}

const ITERATIONS: usize = 1000;
const POPULATION: usize = 500;

// FeasibleCrossoverWrapper global probability parameters
const P_SINGLE_REPLACED: f64 = 0.4;
const P_DOUBLE_FIRST_REPLACED: f64 = 0.6;
const P_DOUBLE_SECOND_REPLACED: f64 = 0.3;
const ARCHIVE_SIZE: usize = 100;
const PARENTS_COUNT: usize = 500;
const G11_EPS: f64 = 0.00015;



@@ 1402,6 1546,56 @@ fn run_nsga_constr<const DIM: usize, const CONSTRAINTS: usize>(
    Ok(())
}

fn run_nsga_improved<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_improved(
        &problem,
        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_improved", &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-Constr handler functions
fn handle_g04_nsga_constr() -> Result<(), Box<dyn std::error::Error>> {
    let problem = problem_g04();


@@ 1491,6 1685,95 @@ fn handle_g24_nsga_constr() -> Result<(), Box<dyn std::error::Error>> {
    run_nsga_constr(problem, config)
}

// NSGA-Improved handler functions
fn handle_g04_nsga_improved() -> 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: 0.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g05_nsga_improved() -> 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: 0.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g06_nsga_improved() -> 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: 0.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g08_nsga_improved() -> 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.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g09_nsga_improved() -> 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: 0.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g11_nsga_improved() -> 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.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g21_nsga_improved() -> 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: 0.1,
    };
    run_nsga_improved(problem, config)
}

fn handle_g24_nsga_improved() -> 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_improved(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();


@@ 1689,11 1972,27 @@ fn main() {

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

    // Run the main logic with increased stack size
    let result = std::thread::Builder::new()
        .stack_size(16 * 1024 * 1024)  // 16MB stack
        .spawn(move || main_with_args(args))
        .expect("Failed to spawn thread")
        .join()
        .expect("Thread panicked");

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

fn main_with_args(args: Vec<String>) -> Result<(), Box<dyn std::error::Error + Send>> {

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



@@ 1722,6 2021,14 @@ fn main() {
        ("nsga_constr", "g11") => handle_g11_nsga_constr(),
        ("nsga_constr", "g21") => handle_g21_nsga_constr(),
        ("nsga_constr", "g24") => handle_g24_nsga_constr(),
        ("nsga_improved", "g04") => handle_g04_nsga_improved(),
        ("nsga_improved", "g05") => handle_g05_nsga_improved(),
        ("nsga_improved", "g06") => handle_g06_nsga_improved(),
        ("nsga_improved", "g08") => handle_g08_nsga_improved(),
        ("nsga_improved", "g09") => handle_g09_nsga_improved(),
        ("nsga_improved", "g11") => handle_g11_nsga_improved(),
        ("nsga_improved", "g21") => handle_g21_nsga_improved(),
        ("nsga_improved", "g24") => handle_g24_nsga_improved(),
        ("nsga_multi", "g04") => handle_nsga_multi_g04(true),
        ("nsga_multi", "g05") => handle_nsga_multi_g05(true),
        ("nsga_multi", "g06") => handle_nsga_multi_g06(true),


@@ 1740,7 2047,7 @@ fn main() {
        ("nsga_multi_noncapped", "g24") => handle_nsga_multi_g24(false),
        (_, _) => {
            eprintln!("Invalid method '{}' or problem '{}'", method, problem);
            eprintln!("Methods: srank, nsga, nsga_constr, nsga_multi, nsga_multi_noncapped");
            eprintln!("Methods: srank, nsga, nsga_constr, nsga_improved, nsga_multi, nsga_multi_noncapped");
            eprintln!("Problems: g04, g05, g06, g08, g09, g11, g21, g24");
            std::process::exit(1);
        }


@@ 1750,4 2057,6 @@ fn main() {
        eprintln!("Error: {}", e);
        std::process::exit(1);
    }

    Ok(())
}

M codes/eoa_lib/src/multi_objective_evolution.rs => codes/eoa_lib/src/multi_objective_evolution.rs +8 -11
@@ 171,7 171,6 @@ impl<'a,
    }
}


#[derive(Debug, Clone, PartialEq)]
pub struct NSGAEvaluation<TOut: PartialEq + Clone, const OBJECTIVES: usize> {
    pub evaluations: [TOut; OBJECTIVES],


@@ 558,9 557,8 @@ pub fn nsga_2<const OBJECTIVES: usize,
    mut evolutionary_strategy: impl FnMut(
        usize,
        &EvolutionStats<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>,
        &EvaluatedPopulation<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>

        // TODO
        &EvaluatedPopulation<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>,
        &mut TCrossover
    ),
    better_than_stats: impl Fn(&TChromosome, &NSGAEvaluation<TResult, OBJECTIVES>, &Option<EvolutionCandidate<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>>) -> bool,
) -> Result<EvolutionResult<TChromosome, [TResult; OBJECTIVES]>, Box<dyn Error>> {


@@ 576,8 574,8 @@ pub fn nsga_2<const OBJECTIVES: usize,
        &MinimizingOperator::new(),
        iterations,
        rng,
        |iteration, stats, population, _, _, _, _, _, _| {
            evolutionary_strategy(iteration, stats, population);
        |iteration, stats, population, _, _, _, crossover, _, _| {
            evolutionary_strategy(iteration, stats, population, crossover);
        },
        better_than_stats
    )?;


@@ 613,9 611,8 @@ pub fn constrained_nsga_2<
    mut evolutionary_strategy: impl FnMut(
        usize,
        &EvolutionStats<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>,
        &EvaluatedPopulation<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>

        // TODO
        &EvaluatedPopulation<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>,
        &mut TCrossover
    ),
    better_than_stats: impl Fn(&TChromosome, &NSGAEvaluation<TResult, OBJECTIVES>, &Option<EvolutionCandidate<TChromosome, NSGAEvaluation<TResult, OBJECTIVES>>>) -> bool,
) -> Result<EvolutionResult<TChromosome, [TResult; OBJECTIVES]>, Box<dyn Error>> {


@@ 631,8 628,8 @@ pub fn constrained_nsga_2<
        &MinimizingOperator::new(),
        iterations,
        rng,
        |iteration, stats, population, _, _, _, _, _, _| {
            evolutionary_strategy(iteration, stats, population);
        |iteration, stats, population, _, _, _, crossover, _, _| {
            evolutionary_strategy(iteration, stats, population, crossover);
        },
        better_than_stats
    )?;