pub mod tsp;
pub mod union_find;
pub mod initializers;
pub mod crossovers;
pub mod binary_string_representation;
pub mod perturbations;
pub mod graph;

use tsp::{NodePermutation, TSPInstance};
use initializers::{MinimumSpanningTreeInitializer, NearestNeighborInitializer, TSPRandomInitializer};
use crossovers::{CycleCrossover, EdgeRecombinationCrossover, NoCrossover, PartiallyMappedCrossover};
use perturbations::{MovePerturbation, Random2OptPerturbation, ReverseSubsequencePerturbation, SwapPerturbation};
use binary_string_representation::TSPBinaryStringWrapper;
use nalgebra::{Dim, Dyn};
use eoa_lib::{
    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;
use std::fs::{File, create_dir_all};
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") {
        Box::new(BufReader::new(GzDecoder::new(file)))
    } else {
        Box::new(BufReader::new(file))
    };

    let mut cities = Vec::new();
    let mut in_coord_section = false;
    let mut dimension = 0;

    for line in reader.lines() {
        let line = line?.trim().to_string();

        if line.starts_with("DIMENSION") {
            dimension = line.split(':').nth(1)
                .or_else(|| line.split_whitespace().nth(1))
                .ok_or("Could not parse dimension")?
                .trim()
                .parse::<usize>()?;
        } else if line == "NODE_COORD_SECTION" {
            in_coord_section = true;
        } else if line == "EOF" {
            break;
        } else if in_coord_section && !line.is_empty() {
            let parts: Vec<&str> = line.split_whitespace().collect();
            if parts.len() >= 3 {
                let x: f64 = parts[1].parse()?;
                let y: f64 = parts[2].parse()?;
                cities.push((x, y));
            }
        }
    }

    if cities.len() != dimension {
        return Err(format!("Expected {} cities, found {}", dimension, cities.len()).into());
    }

    Ok(TSPInstance::new_dyn(cities))
}

#[derive(Debug, Clone)]
struct PlotData {
    best_solution: NodePermutation<Dyn>,
    iterations: Vec<usize>,
    evaluations: Vec<f64>,
    final_cost: f64,
    total_iterations: usize,
    algorithm_name: String,
}

fn extract_evolution_data(
    stats: &EvolutionStats<NodePermutation<Dyn>, f64>,
    final_solution: &NodePermutation<Dyn>,
    final_evaluation: f64,
    final_iteration: usize,
) -> PlotData {
    let mut iterations = Vec::new();
    let mut evaluations = Vec::new();

    for candidate in &stats.best_candidates {
        iterations.push(candidate.evaluation);
        evaluations.push(candidate.evaluated_chromosome.evaluation);
    }

    PlotData {
        best_solution: final_solution.clone(),
        iterations,
        evaluations,
        final_cost: final_evaluation,
        total_iterations: final_iteration,
        algorithm_name: "Evolution Algorithm".to_string(),
    }
}

fn extract_binary_evolution_data(
    stats: &EvolutionStats<BinaryString<Dyn>, f64>,
    final_solution: &NodePermutation<Dyn>,
    final_evaluation: f64,
    final_iteration: usize,
) -> PlotData {
    let mut iterations = Vec::new();
    let mut evaluations = Vec::new();

    for candidate in &stats.best_candidates {
        iterations.push(candidate.evaluation);
        evaluations.push(candidate.evaluated_chromosome.evaluation);
    }

    PlotData {
        best_solution: final_solution.clone(),
        iterations,
        evaluations,
        final_cost: final_evaluation,
        total_iterations: final_iteration,
        algorithm_name: "Evolution Algorithm (Binary)".to_string(),
    }
}

fn extract_local_search_data(
    stats: &LocalSearchStats<NodePermutation<Dyn>, f64>,
    final_solution: &NodePermutation<Dyn>,
    final_evaluation: f64,
    final_cycle: usize,
) -> PlotData {
    let mut iterations = Vec::new();
    let mut evaluations = Vec::new();

    for candidate in stats.candidates() {
        iterations.push(candidate.cycle);
        evaluations.push(candidate.fit);
    }

    // Add final result
    iterations.push(final_cycle);
    evaluations.push(final_evaluation);

    PlotData {
        best_solution: final_solution.clone(),
        iterations,
        evaluations,
        final_cost: final_evaluation,
        total_iterations: final_cycle,
        algorithm_name: "Local Search".to_string(),
    }
}

fn save_results(
    instance: &TSPInstance<Dyn>,
    plot_data: &PlotData,
    base_path: &str,
) -> Result<(), Box<dyn std::error::Error>> {
    // Plot the best solution
    let plot_path = format!("{}.png", base_path);
    instance.draw_solution(&plot_data.best_solution, &plot_path)?;

    // Save statistics to CSV
    let stats_path = format!("{}.csv", base_path);
    let mut stats_file = File::create(&stats_path)?;
    writeln!(stats_file, "fitness_evaluations,evaluation")?;

    for (iteration, evaluation) in plot_data.iterations.iter().zip(plot_data.evaluations.iter()) {
        writeln!(stats_file, "{},{}", iteration, evaluation)?;
    }

    Ok(())
}

fn load_optimal_cost(instance_filename: &str) -> Result<f64, Box<dyn std::error::Error>> {
    let instance_name = std::path::Path::new(instance_filename)
        .file_stem()
        .and_then(|s| s.to_str())
        .ok_or("Could not extract instance name")?
        .trim_end_matches(".tsp");

    let solutions_path = std::path::Path::new(instance_filename)
        .parent()
        .unwrap_or_else(|| std::path::Path::new("."))
        .join("solutions.txt");

    let content = std::fs::read_to_string(solutions_path)?;

    for line in content.lines() {
        let line = line.trim();
        if let Some(colon_pos) = line.find(':') {
            let name = line[..colon_pos].trim();
            if name == instance_name {
                let cost_str = line[colon_pos + 1..].trim();
                return cost_str.parse::<f64>()
                    .map_err(|e| format!("Could not parse cost '{}': {}", cost_str, e).into());
            }
        }
    }

    Err(format!("Optimal cost not found for instance '{}'", instance_name).into())
}

fn run_evolution_algorithm(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 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_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 = 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_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)
}

fn run_evolution_algorithm_binary(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn std::error::Error>> {
    let mut rng = rng();
    let initializer = RandomInitializer::new_binary();
    let output_dimension = instance.dimension();
    let input_dimension = Dyn(output_dimension.value() * (output_dimension.value() - 1) / 2);

    // Create combined perturbation with two mutations wrapped in MutationPerturbation
    let bit_mutation = MutationPerturbation::new(Box::new(BinaryStringBitPerturbation::new(0.1)), 0.2);
    let single_bit_mutation = MutationPerturbation::new(Box::new(BinaryStringSingleBitPerturbation::new()), 0.4);
    let flip1_mutation = MutationPerturbation::new(Box::new(BinaryStringFlipNPerturbation::new(30)), 0.4);
    let flip2_mutation = MutationPerturbation::new(Box::new(BinaryStringFlipNPerturbation::new(20)), 0.4);
    let mut combined_perturbation = CombinedPerturbation::new(vec![
        Box::new(bit_mutation),
        Box::new(single_bit_mutation),
        Box::new(flip1_mutation),
        Box::new(flip2_mutation),
    ]);

    // Set up other components
    let mut crossover = BinaryNPointCrossover::<10, _, _>::new();
    let mut selection = BestSelection::new();
    let mut replacement = TournamentReplacement::new(5, 1.0);
    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(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 = EA_PARENTS_COUNT;
    let result = evolution_algorithm(
        initial_population.clone(),
        parents_count,
        &fitness,
        &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);
            let mut found = false;
            apply_to_perturbations::<_, BinaryStringBitPerturbation<Dyn>>(
                perturbation,
                &mut |p| {
                    found = true;
                    p.p = (0.025 * (1.0 + (iters_since_better as f64 / iters_till_end as f64))).min(0.2);
                }
            );
            assert!(found);

            let mut found = 0;
            MutationPerturbation::apply_to_mutations(
                perturbation,
                &mut |p| {
                    // Do not touch multi bit mutation
                    if found > 0 {
                        p.probability = (0.5 * (1.0 + (iters_since_better as f64 / iters_till_end as f64))).min(1.0);
                    }
                    found += 1;
                }
            );
            assert_eq!(found, 4);
        }
    )?;

    // Extract plotting data
    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,
    );

    Ok(plot_data)
}

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();

    // Create a random initial solution
    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(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_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, 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);
    }

    let instance_name = &args[1];
    let algorithm = &args[2];

    // Load TSP instance
    let filename = format!("instances/{}.tsp.gz", instance_name);
    let instance = load_tsp_instance(&filename)?;
    let dimension = instance.dimension();
    let optimal_cost = load_optimal_cost(&filename)?;

    println!("Loaded {} with {} cities (optimal cost: {})", instance_name, dimension.value(), optimal_cost);

    // Create directory structure: algorithm/instance_name/
    let output_dir = format!("solutions/{}/{}", algorithm, instance_name);
    create_dir_all(&output_dir)?;

    // Generate timestamp for filename
    let now: DateTime<Local> = Local::now();
    let timestamp = now.format("%Y-%m-%d_%H-%M-%S");
    let solution_base_path = format!("{}/{}", output_dir, timestamp);

    // Run the specified algorithm and get plotting data
    let plot_data = match algorithm.as_str() {
        "ea" => {
            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_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_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);
        }
    };

    // Print results
    println!("{} completed in {} iterations", plot_data.algorithm_name, plot_data.total_iterations);
    let initial_cost = plot_data.evaluations[0];
    println!(
        "Initial cost: {:.2} ({:.1}%)",
        initial_cost,
        ((initial_cost - optimal_cost) / optimal_cost) * 100.0
    );
    println!("Final cost: {:.2}", plot_data.final_cost);
    println!("Gap to optimal: {:.2} ({:.1}%)",
        plot_data.final_cost - optimal_cost,
        ((plot_data.final_cost - optimal_cost) / optimal_cost) * 100.0);

    // Save results (plot and CSV)
    save_results(&instance, &plot_data, &solution_base_path)?;

    println!("Created {}.png and {}.csv", solution_base_path, solution_base_path);
    Ok(())
}