~ruther/ctu-fee-eoa

0f5f2903a05e93a8a3e2fd3560f3a33eb1583afc — Rutherther a month ago aaf2bc6
feat: add a lot of algorithm showcases
M .envrc => .envrc +1 -1
@@ 1,3 1,3 @@
use guix
SHELL=$(guix shell -m manifest.scm bash-minimal -- bash -c 'echo $GUIX_ENVIRONMENT')
export LD_LIBRARY_PATH=$SHELL/lib
#export LD_LIBRARY_PATH=$SHELL/lib

M .gitignore => .gitignore +1 -0
@@ 1,2 1,3 @@
*.png
*.svg
target

M codes/Cargo.toml => codes/Cargo.toml +1 -1
@@ 1,3 1,3 @@
[workspace]
resolver = "3"
members = ["eoa_lib", "tsp_hw01"]
members = ["eoa_lib", "tsp_hw01", "tsp_plotter"]

M codes/eoa_lib/src/evolution.rs => codes/eoa_lib/src/evolution.rs +4 -2
@@ 26,7 26,8 @@ pub struct EvolutionResult<TInput, TResult> {
    pub population: EvaluatedPopulation<TInput, TResult>,
    pub stats: EvolutionStats<TInput, TResult>,
    pub best_candidate: EvaluatedChromosome<TInput, TResult>,
    pub iterations: usize
    pub iterations: usize,
    pub evaluations: usize
}

pub fn evolution_algorithm


@@ 163,7 164,8 @@ pub fn evolution_algorithm
        population: current_population,
        best_candidate,
        stats,
        iterations
        iterations,
        evaluations: current_evaluation
    })
}


M codes/eoa_lib/src/perturbation/mod.rs => codes/eoa_lib/src/perturbation/mod.rs +24 -1
@@ 1,7 1,7 @@
use std::{any::Any, borrow::{Borrow, BorrowMut}, marker::PhantomData};

use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, SVector};
use rand::{distr::Distribution, Rng, RngCore};
use rand::{distr::Distribution, Rng, RngCore, prelude::IteratorRandom};
use rand_distr::{uniform, Normal, NormalError, Uniform};

use crate::binary_string::BinaryString;


@@ 432,6 432,29 @@ impl<T: 'static> PerturbationOperator for CombinedPerturbation<T> {
    }
}

pub struct OneOfPerturbation<T> {
    perturbations: Vec<Box<dyn PerturbationOperator<Chromosome = T>>>
}

impl<T> OneOfPerturbation<T> {
    pub fn new(perturbations: Vec<Box<dyn PerturbationOperator<Chromosome = T>>>) -> Self {
        Self {
            perturbations
        }
    }
}

impl<T: 'static> PerturbationOperator for OneOfPerturbation<T> {
    type Chromosome = T;

    fn perturb(&self, chromosome: &mut Self::Chromosome, rng: &mut dyn RngCore) {
        let chosen = (0..self.perturbations.len()).choose(rng);
        if let Some(chosen) = chosen {
            self.perturbations[chosen].perturb(chromosome, rng);
        }
    }
}

#[cfg(test)]
pub mod tests {
    use crate::binary_string::BinaryString;

M codes/tsp_hw01/src/initializers.rs => codes/tsp_hw01/src/initializers.rs +7 -2
@@ 125,7 125,7 @@ where
{
    fn initialize_single(&self, size: D, rng: &mut dyn RngCore) -> NodePermutation<D> {
        let starting_node = rng.random_range(0..size.value());
        self.initialize_from(size, starting_node, rng.random_range(0.0..=0.7), rng)
        self.initialize_from(size, starting_node, 0.0, rng)
    }

    fn initialize(


@@ 170,7 170,12 @@ where

        // 3. Initialize randomly with random p_choose_second (call initialize_single)
        for _ in 0..count {
            population.push(self.initialize_single(size, rng))
            population.push(self.initialize_from(
                size,
                rng.random_range(0..size.value()),
                rng.random_range(0.0..0.5),
                rng
            ));
        }

        population

M codes/tsp_hw01/src/main.rs => codes/tsp_hw01/src/main.rs +585 -43
@@ 13,7 13,7 @@ use perturbations::{MovePerturbation, Random2OptPerturbation, ReverseSubsequence
use binary_string_representation::TSPBinaryStringWrapper;
use nalgebra::{Dim, Dyn};
use eoa_lib::{
    binary_string::BinaryString, comparison::MinimizingOperator, crossover::BinaryNPointCrossover, evolution::{evolution_algorithm, EvolutionStats}, initializer::{Initializer, RandomInitializer}, local_search::{local_search_first_improving, LocalSearchStats}, pairing::AdjacentPairing, perturbation::{apply_to_perturbations, BinaryStringBitPerturbation, BinaryStringFlipNPerturbation, BinaryStringSingleBitPerturbation, CombinedPerturbation, MutationPerturbation, PerturbationOperator}, replacement::{BestReplacement, TournamentReplacement}, selection::{BestSelection, RouletteWheelSelection, TournamentSelection}, terminating::MaximumCyclesTerminatingCondition
    binary_string::BinaryString, comparison::MinimizingOperator, crossover::BinaryNPointCrossover, evolution::{evolution_algorithm, EvolutionStats}, evolutionary_strategy::IdentityStrategy, initializer::{Initializer, RandomInitializer}, local_search::{local_search_first_improving, LocalSearchStats}, pairing::AdjacentPairing, perturbation::{apply_to_perturbations, BinaryStringBitPerturbation, BinaryStringFlipNPerturbation, BinaryStringSingleBitPerturbation, CombinedPerturbation, IdentityPerturbation, MutationPerturbation, OneOfPerturbation, PerturbationOperator}, random_search::random_search, replacement::{BestReplacement, TournamentReplacement}, selection::{BestSelection, RouletteWheelSelection, TournamentSelection}, terminating::MaximumCyclesTerminatingCondition
};
use rand::rng;
use std::env;


@@ 22,6 22,14 @@ use std::io::{BufRead, BufReader, Write};
use flate2::read::GzDecoder;
use chrono::{DateTime, Local};

// Algorithm iteration/cycle constants
const EA_MAX_ITERATIONS: usize = 5000;
const LS_MAX_CYCLES: usize = 250 * 5000 + 500;

// EA population constants
const EA_POPULATION_SIZE: usize = 500;
const EA_PARENTS_COUNT: usize = 250;

fn load_tsp_instance(filename: &str) -> Result<TSPInstance<Dyn>, Box<dyn std::error::Error>> {
    let file = File::open(filename)?;
    let reader: Box<dyn BufRead> = if filename.ends_with(".gz") {


@@ 79,8 87,6 @@ fn extract_evolution_data(
    final_solution: &NodePermutation<Dyn>,
    final_evaluation: f64,
    final_iteration: usize,
    initial_population_size: usize,
    offspring_count: usize,
) -> PlotData {
    let mut iterations = Vec::new();
    let mut evaluations = Vec::new();


@@ 90,11 96,6 @@ fn extract_evolution_data(
        evaluations.push(candidate.evaluated_chromosome.evaluation);
    }

    // Add final result
    let final_fitness_evaluations = initial_population_size + final_iteration * offspring_count;
    iterations.push(final_fitness_evaluations);
    evaluations.push(final_evaluation);

    PlotData {
        best_solution: final_solution.clone(),
        iterations,


@@ 110,8 111,6 @@ fn extract_binary_evolution_data(
    final_solution: &NodePermutation<Dyn>,
    final_evaluation: f64,
    final_iteration: usize,
    initial_population_size: usize,
    offspring_count: usize,
) -> PlotData {
    let mut iterations = Vec::new();
    let mut evaluations = Vec::new();


@@ 121,11 120,6 @@ fn extract_binary_evolution_data(
        evaluations.push(candidate.evaluated_chromosome.evaluation);
    }

    // Add final result
    let final_fitness_evaluations = initial_population_size + final_iteration * offspring_count;
    iterations.push(final_fitness_evaluations);
    evaluations.push(final_evaluation);

    PlotData {
        best_solution: final_solution.clone(),
        iterations,


@@ 231,24 225,151 @@ fn run_evolution_algorithm(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<

    // Set up other components
    let mut crossover = EdgeRecombinationCrossover::new();
    let mut selection = TournamentSelection::new(5, 0.8);
    let mut selection = RouletteWheelSelection::new();
    let mut replacement = BestReplacement::new();
    let mut pairing = AdjacentPairing::new();
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = 500;
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(dimension, population_size, &mut rng);
    // let two_opt = Random2OptPerturbation::new(instance, 10);

    // for individual in initial_population.iter_mut() {
    //     two_opt.perturb(individual, &mut rng);
    // }
    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    // Run evolution algorithm
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,
        instance,
        &mut selection,
        &mut pairing,
        &mut crossover,
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |iteration, stats, _, _, _, _, perturbation, _| {
            let iters_till_end = EA_MAX_ITERATIONS - iteration + 1;
            let iters_since_better =
                iteration - stats.best_candidates.last().map(|c| c.iteration).unwrap_or(0);
            MutationPerturbation::apply_to_mutations(
                perturbation,
                &mut |p| {
                    p.probability = (0.5 * (1.0 + (iters_since_better as f64 / iters_till_end as f64))).min(1.0);
                }
            );
        }
    )?;

    // Extract plotting data
    let plot_data = extract_evolution_data(
        &result.stats,
        &result.best_candidate.chromosome,
        result.best_candidate.evaluation,
        result.iterations,
    );

    Ok(plot_data)
}

fn run_evolution_algorithm_mst(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = MinimumSpanningTreeInitializer::new(instance);
    let dimension = instance.dimension();

    // Create combined perturbation with two mutations wrapped in MutationPerturbation
    let move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.1);
    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.1);
    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.1);
    let mut combined_perturbation = CombinedPerturbation::new(vec![
        Box::new(move_mutation),
        Box::new(swap_mutation),
        Box::new(reverse_mutation),
    ]);

    // Set up other components
    let mut crossover = EdgeRecombinationCrossover::new();
    let mut selection = RouletteWheelSelection::new();
    let mut replacement = BestReplacement::new();
    let mut pairing = AdjacentPairing::new();
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(dimension, population_size, &mut rng);

    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    // Run evolution algorithm
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,
        instance,
        &mut selection,
        &mut pairing,
        &mut crossover,
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |iteration, stats, _, _, _, _, perturbation, _| {
            let iters_till_end = EA_MAX_ITERATIONS - iteration + 1;
            let iters_since_better =
                iteration - stats.best_candidates.last().map(|c| c.iteration).unwrap_or(0);
            MutationPerturbation::apply_to_mutations(
                perturbation,
                &mut |p| {
                    p.probability = (0.5 * (1.0 + (iters_since_better as f64 / iters_till_end as f64))).min(1.0);
                }
            );
        }
    )?;

    // Extract plotting data
    let plot_data = extract_evolution_data(
        &result.stats,
        &result.best_candidate.chromosome,
        result.best_candidate.evaluation,
        result.iterations,
    );

    Ok(plot_data)
}

fn run_evolution_algorithm_nn(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = NearestNeighborInitializer::new(instance);
    let dimension = instance.dimension();

    // Create combined perturbation with two mutations wrapped in MutationPerturbation
    let move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.1);
    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.1);
    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.1);
    let mut combined_perturbation = CombinedPerturbation::new(vec![
        Box::new(move_mutation),
        Box::new(swap_mutation),
        Box::new(reverse_mutation),
    ]);

    // Set up other components
    let mut crossover = EdgeRecombinationCrossover::new();
    let mut selection = RouletteWheelSelection::new();
    let mut replacement = BestReplacement::new();
    let mut pairing = AdjacentPairing::new();
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(dimension, population_size, &mut rng);

    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    // Run evolution algorithm
    let parents_count = 250;
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,


@@ 259,10 380,10 @@ fn run_evolution_algorithm(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        5000, // max iterations
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |iteration, stats, _, _, _, _, perturbation, _| {
            let iters_till_end = 5000 - iteration + 1;
            let iters_till_end = EA_MAX_ITERATIONS - iteration + 1;
            let iters_since_better =
                iteration - stats.best_candidates.last().map(|c| c.iteration).unwrap_or(0);
            MutationPerturbation::apply_to_mutations(


@@ 275,15 396,185 @@ fn run_evolution_algorithm(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<
    )?;

    // Extract plotting data
    let offspring_count = parents_count / 2;
    let initial_population_size = initial_population.iter().count();
    let plot_data = extract_evolution_data(
        &result.stats,
        &result.best_candidate.chromosome,
        result.best_candidate.evaluation,
        result.iterations,
        initial_population_size,
        offspring_count,
    );

    Ok(plot_data)
}

fn run_evolution_algorithm_cx(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Create combined perturbation with lower probabilities and no adaptive changes
    let move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.001);
    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.001);
    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.001);
    let mut combined_perturbation = CombinedPerturbation::new(vec![
        Box::new(move_mutation),
        Box::new(swap_mutation),
        Box::new(reverse_mutation),
    ]);

    // Set up other components with Cycle Crossover
    let mut crossover = CycleCrossover::new();
    let mut selection = RouletteWheelSelection::new();
    let mut replacement = BestReplacement::new();
    let mut pairing = AdjacentPairing::new();
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(dimension, population_size, &mut rng);

    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    // Run evolution algorithm
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,
        instance,
        &mut selection,
        &mut pairing,
        &mut crossover,
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |_iteration, _stats, _, _, _, _, _perturbation, _| {
            // No adaptive perturbation changes - keep probabilities constant
        }
    )?;

    // Extract plotting data
    let plot_data = extract_evolution_data(
        &result.stats,
        &result.best_candidate.chromosome,
        result.best_candidate.evaluation,
        result.iterations,
    );

    Ok(plot_data)
}

fn run_evolution_algorithm_pmx(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Create combined perturbation with lower probabilities and no adaptive changes
    let move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.001);
    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.001);
    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.001);
    let mut combined_perturbation = CombinedPerturbation::new(vec![
        Box::new(move_mutation),
        Box::new(swap_mutation),
        Box::new(reverse_mutation),
    ]);

    // Set up other components with Partially Mapped Crossover
    let mut crossover = PartiallyMappedCrossover::new();
    let mut selection = RouletteWheelSelection::new();
    let mut replacement = BestReplacement::new();
    let mut pairing = AdjacentPairing::new();
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(dimension, population_size, &mut rng);

    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    // Run evolution algorithm
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,
        instance,
        &mut selection,
        &mut pairing,
        &mut crossover,
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |_iteration, _stats, _, _, _, _, _perturbation, _| {
            // No adaptive perturbation changes - keep probabilities constant
        }
    )?;

    // Extract plotting data
    let plot_data = extract_evolution_data(
        &result.stats,
        &result.best_candidate.chromosome,
        result.best_candidate.evaluation,
        result.iterations,
    );

    Ok(plot_data)
}

fn run_evolution_algorithm_erx(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Create combined perturbation with lower probabilities and no adaptive changes
    let move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.001);
    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.001);
    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.001);
    let mut combined_perturbation = CombinedPerturbation::new(vec![
        Box::new(move_mutation),
        Box::new(swap_mutation),
        Box::new(reverse_mutation),
    ]);

    // Set up other components with Edge Recombination Crossover
    let mut crossover = EdgeRecombinationCrossover::new();
    let mut selection = RouletteWheelSelection::new();
    let mut replacement = BestReplacement::new();
    let mut pairing = AdjacentPairing::new();
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(dimension, population_size, &mut rng);

    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    // Run evolution algorithm
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,
        instance,
        &mut selection,
        &mut pairing,
        &mut crossover,
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |_iteration, _stats, _, _, _, _, _perturbation, _| {
            // No adaptive perturbation changes - keep probabilities constant
        }
    )?;

    // Extract plotting data
    let plot_data = extract_evolution_data(
        &result.stats,
        &result.best_candidate.chromosome,
        result.best_candidate.evaluation,
        result.iterations,
    );

    Ok(plot_data)


@@ 315,14 606,14 @@ fn run_evolution_algorithm_binary(instance: &TSPInstance<Dyn>) -> Result<PlotDat
    let better_than_operator = MinimizingOperator::new();

    // Create initial population
    let population_size = 500;
    let population_size = EA_POPULATION_SIZE;
    let initial_population = initializer.initialize(input_dimension, population_size, &mut rng);
    let initial_population = eoa_lib::replacement::Population::from_vec(initial_population);

    let fitness = TSPBinaryStringWrapper::new(instance, input_dimension, output_dimension).unwrap();

    // Run evolution algorithm
    let parents_count = 250;
    let parents_count = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,


@@ 333,10 624,10 @@ fn run_evolution_algorithm_binary(instance: &TSPInstance<Dyn>) -> Result<PlotDat
        &mut combined_perturbation,
        &mut replacement,
        &better_than_operator,
        5000, // max iterations
        EA_MAX_ITERATIONS, // max iterations
        &mut rng,
        |iteration, stats, _, _, _, _, perturbation, _| {
            let iters_till_end = 5000 - iteration + 1;
            let iters_till_end = EA_MAX_ITERATIONS - iteration + 1;
            let iters_since_better =
                iteration - stats.best_candidates.last().map(|c| c.iteration).unwrap_or(0);
            let mut found = false;


@@ 365,22 656,18 @@ fn run_evolution_algorithm_binary(instance: &TSPInstance<Dyn>) -> Result<PlotDat
    )?;

    // Extract plotting data
    let offspring_count = parents_count / 2;
    let initial_population_size = initial_population.iter().count();
    let best_permutation = fitness.to_permutation(&result.best_candidate.chromosome).unwrap();
    let plot_data = extract_binary_evolution_data(
        &result.stats,
        &best_permutation,
        result.best_candidate.evaluation,
        result.iterations,
        initial_population_size,
        offspring_count,
    );

    Ok(plot_data)
}

fn run_local_search(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
fn run_local_search_reverse(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();


@@ 389,8 676,41 @@ fn run_local_search(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std
    let initial_solution = initializer.initialize_single(dimension, &mut rng);

    // Run local search
    let mut perturbation = ReverseSubsequencePerturbation::new();
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = local_search_first_improving(
        instance,
        &mut terminating_condition,
        &mut perturbation,
        &better_than_operator,
        &initial_solution,
        &mut rng,
    )?;

    // Extract plotting data
    let plot_data = extract_local_search_data(
        &result.stats,
        &result.best_candidate.pos,
        result.best_candidate.fit,
        result.cycles,
    );

    Ok(plot_data)
}

fn run_local_search_mst(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = MinimumSpanningTreeInitializer::new(instance);
    let dimension = instance.dimension();

    // Create a MST-based initial solution
    let initial_solution = initializer.initialize_single(dimension, &mut rng);

    // Run local search
    let mut perturbation = MovePerturbation::new();
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(250 * 5000 + 500);
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = local_search_first_improving(


@@ 413,13 733,191 @@ fn run_local_search(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std
    Ok(plot_data)
}

fn run_local_search_nn(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = NearestNeighborInitializer::new(instance);
    let dimension = instance.dimension();

    // Create a nearest neighbor-based initial solution
    let initial_solution = initializer.initialize_single(dimension, &mut rng);

    // Run local search
    let mut perturbation = MovePerturbation::new();
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = local_search_first_improving(
        instance,
        &mut terminating_condition,
        &mut perturbation,
        &better_than_operator,
        &initial_solution,
        &mut rng,
    )?;

    // Extract plotting data
    let plot_data = extract_local_search_data(
        &result.stats,
        &result.best_candidate.pos,
        result.best_candidate.fit,
        result.cycles,
    );

    Ok(plot_data)
}

fn run_local_search_swap(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Create a random initial solution
    let initial_solution = initializer.initialize_single(dimension, &mut rng);

    // Run local search
    let mut perturbation = SwapPerturbation::new();
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = local_search_first_improving(
        instance,
        &mut terminating_condition,
        &mut perturbation,
        &better_than_operator,
        &initial_solution,
        &mut rng,
    )?;

    // Extract plotting data
    let plot_data = extract_local_search_data(
        &result.stats,
        &result.best_candidate.pos,
        result.best_candidate.fit,
        result.cycles,
    );

    Ok(plot_data)
}

fn run_local_search_move(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Create a random initial solution
    let initial_solution = initializer.initialize_single(dimension, &mut rng);

    // Run local search
    let mut perturbation = MovePerturbation::new();
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = local_search_first_improving(
        instance,
        &mut terminating_condition,
        &mut perturbation,
        &better_than_operator,
        &initial_solution,
        &mut rng,
    )?;

    // Extract plotting data
    let plot_data = extract_local_search_data(
        &result.stats,
        &result.best_candidate.pos,
        result.best_candidate.fit,
        result.cycles,
    );

    Ok(plot_data)
}

fn run_local_search_mix(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Create a random initial solution
    let initial_solution = initializer.initialize_single(dimension, &mut rng);

    // Run local search with mixed perturbations
    let mut perturbation = OneOfPerturbation::new(vec![
        Box::new(MovePerturbation::new()),
        Box::new(SwapPerturbation::new()),
        Box::new(ReverseSubsequencePerturbation::new()),
    ]);
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = local_search_first_improving(
        instance,
        &mut terminating_condition,
        &mut perturbation,
        &better_than_operator,
        &initial_solution,
        &mut rng,
    )?;

    // Extract plotting data
    let plot_data = extract_local_search_data(
        &result.stats,
        &result.best_candidate.pos,
        result.best_candidate.fit,
        result.cycles,
    );

    Ok(plot_data)
}

fn run_random_search(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = TSPRandomInitializer::new();
    let dimension = instance.dimension();

    // Run random search
    let mut terminating_condition = MaximumCyclesTerminatingCondition::new(LS_MAX_CYCLES);
    let better_than_operator = MinimizingOperator::new();

    let result = random_search::<
        Dyn,
        NodePermutation<Dyn>,
        f64,
        TSPInstance<Dyn>,
        MaximumCyclesTerminatingCondition,
        IdentityPerturbation<NodePermutation<Dyn>>,
        MinimizingOperator,
        IdentityStrategy,
        TSPRandomInitializer<Dyn>
    >(
        instance,
        &mut terminating_condition,
        &better_than_operator,
        &initializer,
        dimension,
        &mut rng,
    )?;

    // Extract plotting data
    let mut plot_data = extract_local_search_data(
        &result.stats,
        &result.best_candidate.pos,
        result.best_candidate.fit,
        result.cycles,
    );

    // Update algorithm name for random search
    plot_data.algorithm_name = "Random Search".to_string();

    Ok(plot_data)
}

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let args: Vec<String> = env::args().collect();

    if args.len() != 3 {
        eprintln!("Usage: {} <instance_name> <algorithm>", args[0]);
        eprintln!("  instance_name: e.g., kroA100, berlin52, eil51");
        eprintln!("  algorithm: ea (evolution algorithm) or ls (local search)");
        eprintln!("  algorithm: ea, ea_mst, ea_nn, ea_cx, ea_pmx, ea_erx, ls_reverse, ls_swap, ls_move, ls_mix, ls_mst, ls_nn, rs, or ea_binary");
        std::process::exit(1);
    }



@@ 449,16 947,60 @@ fn main() -> Result<(), Box<dyn std::error::Error>> {
            println!("Running Evolution Algorithm...");
            run_evolution_algorithm(&instance)?
        },
        "ea_mst" => {
            println!("Running Evolution Algorithm with MST initialization...");
            run_evolution_algorithm_mst(&instance)?
        },
        "ea_nn" => {
            println!("Running Evolution Algorithm with Nearest Neighbor initialization...");
            run_evolution_algorithm_nn(&instance)?
        },
        "ea_cx" => {
            println!("Running Evolution Algorithm with Cycle Crossover...");
            run_evolution_algorithm_cx(&instance)?
        },
        "ea_pmx" => {
            println!("Running Evolution Algorithm with Partially Mapped Crossover...");
            run_evolution_algorithm_pmx(&instance)?
        },
        "ea_erx" => {
            println!("Running Evolution Algorithm with Edge Recombination Crossover...");
            run_evolution_algorithm_erx(&instance)?
        },
        "ea_binary" => {
            println!("Running Evolution Algorithm (Binary)...");
            run_evolution_algorithm_binary(&instance)?
        },
        "ls" => {
            println!("Running Local Search...");
            run_local_search(&instance)?
        "ls_reverse" => {
            println!("Running Local Search with Reverse Subsequence perturbation...");
            run_local_search_reverse(&instance)?
        },
        "ls_swap" => {
            println!("Running Local Search with Swap perturbation...");
            run_local_search_swap(&instance)?
        },
        "ls_move" => {
            println!("Running Local Search with Move perturbation...");
            run_local_search_move(&instance)?
        },
        "ls_mix" => {
            println!("Running Local Search with mixed perturbations...");
            run_local_search_mix(&instance)?
        },
        "ls_mst" => {
            println!("Running Local Search with MST initialization...");
            run_local_search_mst(&instance)?
        },
        "ls_nn" => {
            println!("Running Local Search with Nearest Neighbor initialization...");
            run_local_search_nn(&instance)?
        },
        "rs" => {
            println!("Running Random Search...");
            run_random_search(&instance)?
        },
        _ => {
            eprintln!("Unknown algorithm: {}. Use 'ea', 'ea_binary', or 'ls'", algorithm);
            eprintln!("Unknown algorithm: {}. Use 'ea', 'ea_mst', 'ea_nn', 'ea_cx', 'ea_pmx', 'ea_erx', 'ls_reverse', 'ls_swap', 'ls_move', 'ls_mix', 'ls_mst', 'ls_nn', 'rs', or 'ea_binary'", algorithm);
            std::process::exit(1);
        }
    };

A codes/tsp_plotter/Cargo.toml => codes/tsp_plotter/Cargo.toml +10 -0
@@ 0,0 1,10 @@
[package]
name = "tsp_plotter"
version = "0.1.0"
edition = "2024"

[dependencies]
plotters = "0.3"
csv = "1.3"
serde = { version = "1.0", features = ["derive"] }
serde_json = "1.0"

A codes/tsp_plotter/src/main.rs => codes/tsp_plotter/src/main.rs +669 -0
@@ 0,0 1,669 @@
use csv::Reader;
use plotters::prelude::*;
use std::collections::HashMap;
use std::fs;
use std::path::Path;
use serde::{Deserialize, Serialize};

#[derive(Debug, Clone, Serialize, Deserialize)]
#[serde(rename_all = "snake_case")]
enum PlotType {
    FitnessEvolution,
    SuccessProbability,
}

#[derive(Debug)]
struct DataPoint {
    evaluations: u32,
    fitness: f64,
    percentage_deviation: f64,
}

#[derive(Debug, Clone, Serialize, Deserialize)]
struct PlotConfig {
    instances: Vec<String>,
    algorithms: Vec<String>,
    group_by_algorithm: bool,
    base_path: String,
    output_path: String,
    optimal_solutions: HashMap<String, f64>,
    targets: Vec<f64>,
    plot_type: PlotType,
    average_targets: bool,
    algorithm_labels: Option<HashMap<String, String>>,
}

impl Default for PlotConfig {
    fn default() -> Self {
        let mut optimal_solutions = HashMap::new();
        optimal_solutions.insert("eil51".to_string(), 426.0);
        optimal_solutions.insert("kroA100".to_string(), 21282.0);
        
        Self {
            instances: vec!["eil51".to_string()],
            algorithms: vec!["ea".to_string(), "ls".to_string(), "rs".to_string()],
            group_by_algorithm: true,
            base_path: "../tsp_hw01/solutions".to_string(),
            output_path: "comparison_eil51.svg".to_string(),
            optimal_solutions,
            targets: vec![1.0, 5.0, 10.0],
            plot_type: PlotType::FitnessEvolution,
            average_targets: false,
            algorithm_labels: None,
        }
    }
}

fn calculate_percentage_deviation(fitness: f64, optimal: f64) -> f64 {
    ((fitness - optimal) / optimal) * 100.0
}

fn get_algorithm_label(algorithm: &str, config: &PlotConfig) -> String {
    config
        .algorithm_labels
        .as_ref()
        .and_then(|labels| labels.get(algorithm))
        .map(|label| label.clone())
        .unwrap_or_else(|| algorithm.to_uppercase())
}

fn create_step_function(data: Vec<DataPoint>) -> Vec<DataPoint> {
    if data.is_empty() {
        return data;
    }
    
    let mut result = Vec::new();
    
    for i in 0..data.len() {
        let current_point = &data[i];
        
        // Add the actual data point (the vertical part of the step)
        result.push(DataPoint {
            evaluations: current_point.evaluations,
            fitness: current_point.fitness,
            percentage_deviation: current_point.percentage_deviation,
        });
        
        // If this is not the last point, add a horizontal step to the next evaluation
        if i + 1 < data.len() {
            let next_evaluation = data[i + 1].evaluations;
            // Add a point just before the next evaluation with current fitness value
            // This creates the horizontal step line
            result.push(DataPoint {
                evaluations: next_evaluation,
                fitness: current_point.fitness,
                percentage_deviation: current_point.percentage_deviation,
            });
        }
    }
    
    result
}

#[derive(Debug)]
struct ProbabilityPoint {
    evaluations: u32,
    probability: f64,
}

fn calculate_success_probability(
    algorithm_data: &HashMap<String, Vec<DataPoint>>,
    target_percentage: f64,
) -> Vec<ProbabilityPoint> {
    if algorithm_data.is_empty() {
        return Vec::new();
    }
    
    // Collect all unique evaluation points across all runs
    let mut all_evaluations = std::collections::BTreeSet::new();
    for (_, points) in algorithm_data {
        for point in points {
            all_evaluations.insert(point.evaluations);
        }
    }
    
    let mut probability_points = Vec::new();
    
    for &evaluation in &all_evaluations {
        let total_runs = algorithm_data.len();
        let mut successful_runs = 0;
        
        // For each run, check if it has achieved the target at this evaluation
        for (_, points) in algorithm_data {
            // Find the best performance achieved up to this evaluation
            let best_percentage = points
                .iter()
                .filter(|p| p.evaluations <= evaluation)
                .map(|p| p.percentage_deviation)
                .fold(f64::INFINITY, f64::min);
            
            if best_percentage <= target_percentage {
                successful_runs += 1;
            }
        }
        
        let probability = successful_runs as f64 / total_runs as f64;
        probability_points.push(ProbabilityPoint {
            evaluations: evaluation,
            probability,
        });
    }
    
    probability_points
}

fn calculate_averaged_success_probability(
    algorithm_data: &HashMap<String, Vec<DataPoint>>,
    targets: &[f64],
) -> Vec<ProbabilityPoint> {
    if algorithm_data.is_empty() || targets.is_empty() {
        return Vec::new();
    }
    
    // Calculate probability for each target
    let target_probabilities: Vec<Vec<ProbabilityPoint>> = targets
        .iter()
        .map(|&target| calculate_success_probability(algorithm_data, target))
        .collect();
    
    if target_probabilities.is_empty() {
        return Vec::new();
    }
    
    // Collect all unique evaluation points across all targets
    let mut all_evaluations = std::collections::BTreeSet::new();
    for prob_data in &target_probabilities {
        for point in prob_data {
            all_evaluations.insert(point.evaluations);
        }
    }
    
    let mut averaged_points = Vec::new();
    
    for &evaluation in &all_evaluations {
        let mut total_probability = 0.0;
        let mut count = 0;
        
        // Average probabilities across all targets for this evaluation
        for prob_data in &target_probabilities {
            if let Some(point) = prob_data.iter().find(|p| p.evaluations == evaluation) {
                total_probability += point.probability;
                count += 1;
            }
        }
        
        if count > 0 {
            let averaged_probability = total_probability / count as f64;
            averaged_points.push(ProbabilityPoint {
                evaluations: evaluation,
                probability: averaged_probability,
            });
        }
    }
    
    averaged_points
}

fn read_csv_file(file_path: &Path, optimal_solution: f64) -> Result<Vec<DataPoint>, Box<dyn std::error::Error>> {
    let mut reader = Reader::from_path(file_path)?;
    let mut data = Vec::new();

    for result in reader.records() {
        let record = result?;
        if record.len() >= 2 {
            let evaluations: u32 = record[0].parse()?;
            let fitness: f64 = record[1].parse()?;
            let percentage_deviation = calculate_percentage_deviation(fitness, optimal_solution);
            data.push(DataPoint { evaluations, fitness, percentage_deviation });
        }
    }

    // Add a point at evaluation 1 with the same fitness as the first point
    if !data.is_empty() {
        let first_fitness = data[0].fitness;
        let first_percentage = data[0].percentage_deviation;
        data.insert(0, DataPoint { 
            evaluations: 1, 
            fitness: first_fitness,
            percentage_deviation: first_percentage
        });
    }

    Ok(data)
}

fn read_all_csv_files(directory: &Path, optimal_solution: f64) -> Result<HashMap<String, Vec<DataPoint>>, Box<dyn std::error::Error>> {
    let mut all_data = HashMap::new();

    if directory.is_dir() {
        for entry in fs::read_dir(directory)? {
            let entry = entry?;
            let path = entry.path();

            if path.extension().and_then(|s| s.to_str()) == Some("csv") {
                let file_name = path.file_stem()
                    .and_then(|s| s.to_str())
                    .unwrap_or("unknown")
                    .to_string();

                match read_csv_file(&path, optimal_solution) {
                    Ok(data) => {
                        all_data.insert(file_name, data);
                    }
                    Err(e) => {
                        eprintln!("Error reading {}: {}", path.display(), e);
                    }
                }
            }
        }
    }

    Ok(all_data)
}

#[derive(Debug)]
struct PlotData {
    data: HashMap<String, HashMap<String, Vec<DataPoint>>>,
}

fn read_plot_data(config: &PlotConfig) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut data = HashMap::new();
    let base_path = Path::new(&config.base_path);

    for algorithm in &config.algorithms {
        let mut algorithm_data = HashMap::new();
        
        for instance in &config.instances {
            let path = base_path.join(format!("{}/{}", algorithm, instance));
            
            if path.exists() {
                let optimal_solution = config.optimal_solutions.get(instance)
                    .copied()
                    .unwrap_or_else(|| {
                        eprintln!("Warning: No optimal solution found for instance '{}', using default value 1000.0", instance);
                        1000.0
                    });
                
                let instance_data = read_all_csv_files(&path, optimal_solution)?;
                
                // Apply step function to create proper step visualization
                let mut step_data = HashMap::new();
                for (run_key, points) in instance_data {
                    let step_points = create_step_function(points);
                    step_data.insert(run_key, step_points);
                }
                
                algorithm_data.extend(step_data);
            }
        }
        
        if !algorithm_data.is_empty() {
            data.insert(algorithm.clone(), algorithm_data);
        }
    }

    Ok(PlotData { data })
}

fn get_color_palette() -> Vec<RGBColor> {
    vec![
        BLUE, RED, GREEN, CYAN, MAGENTA, 
        RGBColor(255, 165, 0), // Orange
        RGBColor(128, 0, 128),  // Purple
        RGBColor(255, 192, 203), // Pink
        RGBColor(165, 42, 42),   // Brown
        RGBColor(128, 128, 128), // Gray
    ]
}

fn create_plot(plot_data: &PlotData, config: &PlotConfig) -> Result<(), Box<dyn std::error::Error>> {
    // Create plots directory if it doesn't exist
    fs::create_dir_all("plots")?;
    
    let output_path = Path::new("plots").join(&config.output_path);
    
    let root = SVGBackend::new(&output_path, (1024, 768)).into_drawing_area();
    root.fill(&WHITE)?;

    let mut min_evaluations = u32::MAX;
    let mut max_evaluations = 0u32;

    let mut min_percentage = f64::MAX;
    let mut max_percentage = f64::MIN;

    // Find min/max across all data
    for (_, algorithm_data) in &plot_data.data {
        for (_, points) in algorithm_data {
            for point in points {
                min_evaluations = min_evaluations.min(point.evaluations);
                max_evaluations = max_evaluations.max(point.evaluations);
                min_percentage = min_percentage.min(point.percentage_deviation);
                max_percentage = max_percentage.max(point.percentage_deviation);
            }
        }
    }

    println!("Percentage deviation range: {:.2}% to {:.2}%", min_percentage, max_percentage);

    // Use actual data range with some padding for better visualization
    let padding_factor = 0.1;
    let range = max_percentage - min_percentage;
    let y_min = (min_percentage - range * padding_factor).max(0.1); // Ensure y_min > 0 for log scale
    let y_max = max_percentage + range * padding_factor;

    let mut chart = ChartBuilder::on(&root)
        .margin(10)
        .x_label_area_size(50)
        .y_label_area_size(90)
        .build_cartesian_2d(
            (min_evaluations as f64..max_evaluations as f64).log_scale(),
            (y_min..y_max).log_scale(),
        )?;

    chart
        .configure_mesh()
        .x_desc("Evaluations")
        .y_desc("Percentage deviation from optimal (%)")
        .y_labels(10)
        .x_label_formatter(&|x| {
            let power = x.log10().round() as i32;
            match power {
                0 => "10⁰".to_string(),
                1 => "10¹".to_string(),
                2 => "10²".to_string(),
                3 => "10³".to_string(),
                4 => "10⁴".to_string(),
                5 => "10⁵".to_string(),
                6 => "10⁶".to_string(),
                7 => "10⁷".to_string(),
                8 => "10⁸".to_string(),
                9 => "10⁹".to_string(),
                _ => format!("10^{}", power),
            }
        })
        .y_label_formatter(&|y| {
            if *y >= 100.0 {
                format!("{:.0}%", *y)
            } else if *y >= 10.0 {
                format!("{:.1}%", *y)
            } else {
                format!("{:.2}%", *y)
            }
        })
        .x_max_light_lines(8)
        .y_max_light_lines(12)
        .axis_desc_style(("sans-serif", 18))
        .label_style(("sans-serif", 16))
        .draw()?;

    let colors = get_color_palette();
    let mut color_index = 0;
    let mut legend_added = HashMap::new();

    if config.group_by_algorithm {
        // Group by algorithm: EA, LS, RS each get one color
        for algorithm in &config.algorithms {
            if let Some(algorithm_data) = plot_data.data.get(algorithm) {
                let color = colors[color_index % colors.len()];
                color_index += 1;

                for (_, points) in algorithm_data {
                    let series = chart
                        .draw_series(LineSeries::new(
                            points.iter().map(|p| (p.evaluations as f64, p.percentage_deviation)),
                            &color,
                        ))?;

                    if !legend_added.contains_key(algorithm) {
                        let label = get_algorithm_label(algorithm, config);
                        series.label(label).legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
                        legend_added.insert(algorithm.clone(), true);
                    }
                }
            }
        }
    } else {
        // Group by instance: Each algorithm-instance combination gets different color
        for algorithm in &config.algorithms {
            if let Some(algorithm_data) = plot_data.data.get(algorithm) {
                for instance in &config.instances {
                    for (run_key, points) in algorithm_data {
                        if run_key.contains(instance) || config.instances.len() == 1 {
                            let color = colors[color_index % colors.len()];
                            color_index += 1;

                            let series = chart
                                .draw_series(LineSeries::new(
                                    points.iter().map(|p| (p.evaluations as f64, p.percentage_deviation)),
                                    &color,
                                ))?;

                            let algo_label = get_algorithm_label(algorithm, config);
                            let label = format!("{} {}", algo_label, instance);
                            if !legend_added.contains_key(&label) {
                                series.label(&label).legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
                                legend_added.insert(label, true);
                            }
                        }
                    }
                }
            }
        }
    }

    // Draw horizontal target lines
    for &target in &config.targets {
        if target >= y_min && target <= y_max {
            let series = chart
                .draw_series(LineSeries::new(
                    vec![(min_evaluations as f64, target), (max_evaluations as f64, target)],
                    BLACK.stroke_width(2),
                ))?;
                
            series
                .label(&format!("{}% target", target))
                .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], BLACK));
        }
    }

    chart.configure_series_labels()
        .background_style(&WHITE)
        .border_style(&BLACK)
        .position(SeriesLabelPosition::LowerLeft)
        .draw()?;
    root.present()?;

    println!("Plot saved to: {}", output_path.display());
    Ok(())
}

fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) -> Result<(), Box<dyn std::error::Error>> {
    // Create plots directory if it doesn't exist
    fs::create_dir_all("plots")?;
    
    let output_path = Path::new("plots").join(&config.output_path);
    
    let root = SVGBackend::new(&output_path, (1024, 768)).into_drawing_area();
    root.fill(&WHITE)?;

    let mut min_evaluations = u32::MAX;
    let mut max_evaluations = 0u32;

    // Find min/max evaluations across all data
    for (_, algorithm_data) in &plot_data.data {
        for (_, points) in algorithm_data {
            for point in points {
                min_evaluations = min_evaluations.min(point.evaluations);
                max_evaluations = max_evaluations.max(point.evaluations);
            }
        }
    }

    let mut chart = ChartBuilder::on(&root)
        .margin(10)
        .x_label_area_size(50)
        .y_label_area_size(90)
        .build_cartesian_2d(
            (min_evaluations as f64..max_evaluations as f64).log_scale(),
            0.0f64..1.0f64,
        )?;

    chart
        .configure_mesh()
        .x_desc("Evaluations")
        .y_desc("Probability of Success")
        .y_labels(11)
        .x_label_formatter(&|x| {
            let power = x.log10().round() as i32;
            match power {
                0 => "10⁰".to_string(),
                1 => "10¹".to_string(),
                2 => "10²".to_string(),
                3 => "10³".to_string(),
                4 => "10⁴".to_string(),
                5 => "10⁵".to_string(),
                6 => "10⁶".to_string(),
                7 => "10⁷".to_string(),
                8 => "10⁸".to_string(),
                9 => "10⁹".to_string(),
                _ => format!("10^{}", power),
            }
        })
        .y_label_formatter(&|y| {
            format!("{:.1}", *y)
        })
        .x_max_light_lines(8)
        .y_max_light_lines(0)
        .axis_desc_style(("sans-serif", 18))
        .label_style(("sans-serif", 16))
        .draw()?;

    let colors = get_color_palette();
    let mut color_index = 0;

    // Plot probability curves
    if config.average_targets {
        // Plot one averaged curve per algorithm
        for algorithm in &config.algorithms {
            if let Some(algorithm_data) = plot_data.data.get(algorithm) {
                let probability_data = calculate_averaged_success_probability(algorithm_data, &config.targets);
                
                if !probability_data.is_empty() {
                    let color = colors[color_index % colors.len()];
                    color_index += 1;

                    let series = chart
                        .draw_series(LineSeries::new(
                            probability_data.iter().map(|p| (p.evaluations as f64, p.probability)),
                            &color,
                        ))?;

                    let algo_label = get_algorithm_label(algorithm, config);
                    let label = format!("{} (avg)", algo_label);
                    series
                        .label(&label)
                        .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
                }
            }
        }
    } else {
        // Plot individual curves for each algorithm and target combination
        for algorithm in &config.algorithms {
            if let Some(algorithm_data) = plot_data.data.get(algorithm) {
                for &target in &config.targets {
                    let probability_data = calculate_success_probability(algorithm_data, target);
                    
                    if !probability_data.is_empty() {
                        let color = colors[color_index % colors.len()];
                        color_index += 1;

                        let series = chart
                            .draw_series(LineSeries::new(
                                probability_data.iter().map(|p| (p.evaluations as f64, p.probability)),
                                &color,
                            ))?;

                        let algo_label = get_algorithm_label(algorithm, config);
                        let label = format!("{} - {}%", algo_label, target);
                        series
                            .label(&label)
                            .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
                    }
                }
            }
        }
    }

    chart.configure_series_labels().background_style(&WHITE).border_style(&BLACK).draw()?;
    root.present()?;

    println!("Probability plot saved to: {}", output_path.display());
    Ok(())
}

fn load_config(config_path: &str) -> Result<PlotConfig, Box<dyn std::error::Error>> {
    
    if Path::new(config_path).exists() {
        let config_str = fs::read_to_string(config_path)?;
        let config: PlotConfig = serde_json::from_str(&config_str)?;
        println!("Loaded configuration from {}", config_path);
        Ok(config)
    } else {
        let config = PlotConfig::default();
        
        // Create default config file
        let config_str = serde_json::to_string_pretty(&config)?;
        fs::write(config_path, config_str)?;
        println!("Created default configuration file: {}", config_path);
        
        Ok(config)
    }
}

fn main() -> Result<(), Box<dyn std::error::Error>> {
    let args: Vec<String> = std::env::args().collect();
    
    if args.len() != 2 {
        eprintln!("Usage: {} <config_file.json>", args[0]);
        eprintln!("Example: {} plot_config.json", args[0]);
        std::process::exit(1);
    }
    
    let config_path = &args[1];
    let config = load_config(config_path)?;
    
    println!("Configuration:");
    println!("  Instances: {:?}", config.instances);
    println!("  Algorithms: {:?}", config.algorithms);
    println!("  Group by algorithm: {}", config.group_by_algorithm);
    println!("  Base path: {}", config.base_path);
    println!("  Output path: {}", config.output_path);
    println!("  Plot type: {:?}", config.plot_type);
    println!("  Average targets: {}", config.average_targets);
    println!();

    let base_directory = Path::new(&config.base_path);
    println!("Reading CSV files from: {}", base_directory.display());

    let plot_data = read_plot_data(&config)?;

    if plot_data.data.is_empty() {
        eprintln!("No CSV files found for any configured algorithms and instances");
        return Ok(());
    }

    for (algorithm, algorithm_data) in &plot_data.data {
        println!("Found {} CSV files for algorithm: {}", algorithm_data.len(), algorithm);
    }

    match config.plot_type {
        PlotType::FitnessEvolution => {
            create_plot(&plot_data, &config)?;
        }
        PlotType::SuccessProbability => {
            create_success_probability_plot(&plot_data, &config)?;
        }
    }

    Ok(())
}