From 8315b8f68fee89eef358a5209cc43b03855e64bd Mon Sep 17 00:00:00 2001
From: Rutherther <rutherther@ditigal.xyz>
Date: Sun, 2 Nov 2025 22:32:46 +0100
Subject: [PATCH] fix(tsp): few more tweaks

---
 codes/tsp_hw01/src/main.rs    |  30 +--
 codes/tsp_plotter/src/main.rs | 342 +++++++++++++++++++++++++++++-----
 2 files changed, 315 insertions(+), 57 deletions(-)

diff --git a/codes/tsp_hw01/src/main.rs b/codes/tsp_hw01/src/main.rs
index 4068039322187e8a8b90c045aa529a67d3cf87b8..61b524cc69ad0d4cf2250484d83cb1900694420e 100644
--- a/codes/tsp_hw01/src/main.rs
+++ b/codes/tsp_hw01/src/main.rs
@@ -440,9 +440,9 @@ fn run_evolution_algorithm_cx(instance: &TSPInstance<Dyn>) -> Result<PlotData, B
     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 move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.000);
+    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.000);
+    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.000);
     let mut combined_perturbation = CombinedPerturbation::new(vec![
         Box::new(move_mutation),
         Box::new(swap_mutation),
@@ -498,9 +498,9 @@ fn run_evolution_algorithm_pmx(instance: &TSPInstance<Dyn>) -> Result<PlotData,
     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 move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.000);
+    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.000);
+    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.000);
     let mut combined_perturbation = CombinedPerturbation::new(vec![
         Box::new(move_mutation),
         Box::new(swap_mutation),
@@ -556,9 +556,9 @@ fn run_evolution_algorithm_erx(instance: &TSPInstance<Dyn>) -> Result<PlotData,
     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 move_mutation = MutationPerturbation::new(Box::new(MovePerturbation::new()), 0.000);
+    let swap_mutation = MutationPerturbation::new(Box::new(SwapPerturbation::new()), 0.000);
+    let reverse_mutation = MutationPerturbation::new(Box::new(ReverseSubsequencePerturbation::new()), 0.000);
     let mut combined_perturbation = CombinedPerturbation::new(vec![
         Box::new(move_mutation),
         Box::new(swap_mutation),
@@ -737,7 +737,11 @@ fn run_local_search_mst(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn
     let initial_solution = initializer.initialize_single(dimension, &mut rng);
 
     // Run local search
-    let mut perturbation = MovePerturbation::new();
+    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();
 
@@ -770,7 +774,11 @@ fn run_local_search_nn(instance: &TSPInstance<Dyn>) -> Result<PlotData, Box<dyn
     let initial_solution = initializer.initialize_single(dimension, &mut rng);
 
     // Run local search
-    let mut perturbation = MovePerturbation::new();
+    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();
 
diff --git a/codes/tsp_plotter/src/main.rs b/codes/tsp_plotter/src/main.rs
index 0885cc1133c8a91cdcb63aad29644c73e3cd231c..a213c4f05c3dcf7e2175b151ce3a5294dbc39792 100644
--- a/codes/tsp_plotter/src/main.rs
+++ b/codes/tsp_plotter/src/main.rs
@@ -13,7 +13,7 @@ enum PlotType {
     SuccessProbability,
 }
 
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 struct DataPoint {
     evaluations: u32,
     fitness: f64,
@@ -31,6 +31,18 @@ struct PlotConfig {
     plot_type: PlotType,
     average_targets: bool,
     algorithm_labels: Option<HashMap<String, String>>,
+    #[serde(default = "default_show_std_dev")]
+    show_std_dev: bool,
+    #[serde(default = "default_average_runs")]
+    average_runs: bool,
+}
+
+fn default_show_std_dev() -> bool {
+    true
+}
+
+fn default_average_runs() -> bool {
+    false
 }
 
 impl Default for PlotConfig {
@@ -45,6 +57,8 @@ impl Default for PlotConfig {
             plot_type: PlotType::FitnessEvolution,
             average_targets: false,
             algorithm_labels: None,
+            show_std_dev: true,
+            average_runs: false,
         }
     }
 }
@@ -127,13 +141,13 @@ fn create_step_function(data: Vec<DataPoint>) -> Vec<DataPoint> {
     result
 }
 
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 struct ProbabilityPoint {
     evaluations: u32,
     probability: f64,
 }
 
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 struct ProbabilityPointWithDeviation {
     evaluations: u32,
     probability: f64,
@@ -142,6 +156,16 @@ struct ProbabilityPointWithDeviation {
     upper_bound: f64,
 }
 
+#[derive(Debug, Clone)]
+struct DataPointWithDeviation {
+    evaluations: u32,
+    fitness: f64,
+    percentage_deviation: f64,
+    std_dev: f64,
+    lower_bound: f64,
+    upper_bound: f64,
+}
+
 fn calculate_success_probability(
     algorithm_data: &HashMap<String, Vec<DataPoint>>,
     target_percentage: f64,
@@ -150,9 +174,35 @@ fn calculate_success_probability(
         return Vec::new();
     }
 
+    // Find the maximum evaluation point across all data
+    let max_evaluation = algorithm_data
+        .values()
+        .flat_map(|points| points.iter())
+        .map(|p| p.evaluations)
+        .max()
+        .unwrap_or(0);
+
+    // Extend all data to the maximum evaluation point
+    let mut extended_algorithm_data = HashMap::new();
+    for (run_key, points) in algorithm_data {
+        let mut extended_points = points.clone();
+        if !extended_points.is_empty() {
+            extended_points.sort_by_key(|point| point.evaluations);
+            let last_point = extended_points.last().unwrap();
+            if last_point.evaluations < max_evaluation {
+                extended_points.push(DataPoint {
+                    evaluations: max_evaluation,
+                    fitness: last_point.fitness,
+                    percentage_deviation: last_point.percentage_deviation,
+                });
+            }
+        }
+        extended_algorithm_data.insert(run_key.clone(), extended_points);
+    }
+
     // Collect all unique evaluation points across all runs
     let mut all_evaluations = std::collections::BTreeSet::new();
-    for (_, points) in algorithm_data {
+    for (_, points) in &extended_algorithm_data {
         for point in points {
             all_evaluations.insert(point.evaluations);
         }
@@ -161,11 +211,11 @@ fn calculate_success_probability(
     let mut probability_points = Vec::new();
 
     for &evaluation in &all_evaluations {
-        let total_runs = algorithm_data.len();
+        let total_runs = extended_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 {
+        for (_, points) in &extended_algorithm_data {
             // Find the best performance achieved up to this evaluation
             let best_percentage = points
                 .iter()
@@ -252,6 +302,107 @@ fn calculate_averaged_success_probability_with_deviation(
     averaged_points
 }
 
+fn calculate_averaged_fitness_with_deviation(
+    algorithm_data: &HashMap<String, Vec<DataPoint>>,
+) -> Vec<DataPointWithDeviation> {
+    if algorithm_data.is_empty() {
+        return Vec::new();
+    }
+
+    // Find the maximum evaluation point across all data
+    let max_evaluation = algorithm_data
+        .values()
+        .flat_map(|points| points.iter())
+        .map(|p| p.evaluations)
+        .max()
+        .unwrap_or(0);
+
+    // Extend all data to the maximum evaluation point
+    let mut extended_algorithm_data = HashMap::new();
+    for (run_key, points) in algorithm_data {
+        let mut extended_points = points.clone();
+        if !extended_points.is_empty() {
+            extended_points.sort_by_key(|point| point.evaluations);
+            let last_point = extended_points.last().unwrap();
+            if last_point.evaluations < max_evaluation {
+                extended_points.push(DataPoint {
+                    evaluations: max_evaluation,
+                    fitness: last_point.fitness,
+                    percentage_deviation: last_point.percentage_deviation,
+                });
+            }
+        }
+        extended_algorithm_data.insert(run_key.clone(), extended_points);
+    }
+
+    // Collect all unique evaluation points across all runs
+    let mut all_evaluations = std::collections::BTreeSet::new();
+    for (_, points) in &extended_algorithm_data {
+        for point in points {
+            all_evaluations.insert(point.evaluations);
+        }
+    }
+
+    let mut averaged_points = Vec::new();
+
+    for &evaluation in &all_evaluations {
+        let mut fitness_values = Vec::new();
+        let mut percentage_values = Vec::new();
+
+        // Collect fitness and percentage values at this evaluation point from all runs
+        for (_, points) in &extended_algorithm_data {
+            // Find the best (minimum) percentage deviation 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);
+            
+            // Find the corresponding fitness value
+            let best_fitness = points
+                .iter()
+                .filter(|p| p.evaluations <= evaluation)
+                .min_by(|a, b| a.percentage_deviation.partial_cmp(&b.percentage_deviation).unwrap())
+                .map(|p| p.fitness)
+                .unwrap_or(f64::INFINITY);
+
+            if best_percentage != f64::INFINITY {
+                fitness_values.push(best_fitness);
+                percentage_values.push(best_percentage);
+            }
+        }
+
+        if !fitness_values.is_empty() {
+            // Calculate means
+            let mean_fitness = fitness_values.iter().sum::<f64>() / fitness_values.len() as f64;
+            let mean_percentage = percentage_values.iter().sum::<f64>() / percentage_values.len() as f64;
+            
+            // Calculate standard deviation for percentage deviation
+            let variance = percentage_values.iter()
+                .map(|p| (p - mean_percentage).powi(2))
+                .sum::<f64>() / percentage_values.len() as f64;
+            let std_dev = variance.sqrt();
+            
+            // Calculate bounds (mean ± 1 standard deviation, clamped to reasonable values)
+            let lower_bound = (mean_percentage - std_dev).max(0.1); // Keep above 0.1% for log scale
+            let upper_bound = mean_percentage + std_dev;
+            
+            averaged_points.push(DataPointWithDeviation {
+                evaluations: evaluation,
+                fitness: mean_fitness,
+                percentage_deviation: mean_percentage,
+                std_dev,
+                lower_bound,
+                upper_bound,
+            });
+        }
+    }
+
+    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)?;
@@ -451,17 +602,90 @@ fn create_plot(plot_data: &PlotData, config: &PlotConfig) -> Result<(), Box<dyn
                 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 config.average_runs {
+                    // Calculate averaged fitness data with deviation
+                    let averaged_data = calculate_averaged_fitness_with_deviation(algorithm_data);
+                    
+                    if !averaged_data.is_empty() {
+                        // Conditionally draw standard deviation bands
+                        if config.show_std_dev {
+                            // Create transparent confidence band
+                            let transparent_color = color.mix(0.3); // 30% opacity
+                            
+                            // Create upper and lower bound points for the filled area
+                            let upper_points: Vec<(f64, f64)> = averaged_data
+                                .iter()
+                                .map(|p| (p.evaluations as f64, p.upper_bound))
+                                .collect();
+                            
+                            let mut lower_points: Vec<(f64, f64)> = averaged_data
+                                .iter()
+                                .map(|p| (p.evaluations as f64, p.lower_bound))
+                                .collect();
+                            
+                            // Reverse the lower points to create a closed polygon
+                            lower_points.reverse();
+                            
+                            // Combine upper and lower points to form a polygon
+                            let mut polygon_points = upper_points;
+                            polygon_points.extend(lower_points);
+                            
+                            // Draw the filled confidence band
+                            if polygon_points.len() > 2 {
+                                chart.draw_series(std::iter::once(Polygon::new(
+                                    polygon_points,
+                                    transparent_color.filled(),
+                                )))?;
+                            }
+                        }
 
-                    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);
+                        // Draw the main averaged line on top of the confidence band
+                        let series = chart
+                            .draw_series(LineSeries::new(
+                                averaged_data.iter().map(|p| (p.evaluations as f64, p.percentage_deviation)),
+                                &color,
+                            ))?;
+
+                        if !legend_added.contains_key(algorithm) {
+                            let label = get_algorithm_label(algorithm, config);
+                            let final_label = if config.show_std_dev {
+                                format!("{} (avg ± σ)", label)
+                            } else {
+                                format!("{} (avg)", label)
+                            };
+                            series.label(final_label).legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
+                            legend_added.insert(algorithm.clone(), true);
+                        }
+                    }
+                } else {
+                    // Original individual run plotting
+                    for (_, points) in algorithm_data {
+                        let mut extended_points = points.clone();
+                        
+                        // Extend data to max evaluation inline
+                        if !extended_points.is_empty() {
+                            extended_points.sort_by_key(|point| point.evaluations);
+                            let last_point = extended_points.last().unwrap();
+                            if last_point.evaluations < max_evaluations {
+                                extended_points.push(DataPoint {
+                                    evaluations: max_evaluations,
+                                    fitness: last_point.fitness,
+                                    percentage_deviation: last_point.percentage_deviation,
+                                });
+                            }
+                        }
+                        
+                        let series = chart
+                            .draw_series(LineSeries::new(
+                                extended_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);
+                        }
                     }
                 }
             }
@@ -476,9 +700,24 @@ fn create_plot(plot_data: &PlotData, config: &PlotConfig) -> Result<(), Box<dyn
                             let color = colors[color_index % colors.len()];
                             color_index += 1;
 
+                            let mut extended_points = points.clone();
+                            
+                            // Extend data to max evaluation inline
+                            if !extended_points.is_empty() {
+                                extended_points.sort_by_key(|point| point.evaluations);
+                                let last_point = extended_points.last().unwrap();
+                                if last_point.evaluations < max_evaluations {
+                                    extended_points.push(DataPoint {
+                                        evaluations: max_evaluations,
+                                        fitness: last_point.fitness,
+                                        percentage_deviation: last_point.percentage_deviation,
+                                    });
+                                }
+                            }
+                            
                             let series = chart
                                 .draw_series(LineSeries::new(
-                                    points.iter().map(|p| (p.evaluations as f64, p.percentage_deviation)),
+                                    extended_points.iter().map(|p| (p.evaluations as f64, p.percentage_deviation)),
                                     &color,
                                 ))?;
 
@@ -596,33 +835,36 @@ fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) ->
                     let color = colors[color_index % colors.len()];
                     color_index += 1;
 
-                    // Create transparent confidence band
-                    let transparent_color = color.mix(0.3); // 30% opacity
-                    
-                    // Create upper and lower bound points for the filled area
-                    let upper_points: Vec<(f64, f64)> = probability_data
-                        .iter()
-                        .map(|p| (p.evaluations as f64, p.upper_bound))
-                        .collect();
-                    
-                    let mut lower_points: Vec<(f64, f64)> = probability_data
-                        .iter()
-                        .map(|p| (p.evaluations as f64, p.lower_bound))
-                        .collect();
-                    
-                    // Reverse the lower points to create a closed polygon
-                    lower_points.reverse();
-                    
-                    // Combine upper and lower points to form a polygon
-                    let mut polygon_points = upper_points;
-                    polygon_points.extend(lower_points);
-                    
-                    // Draw the filled confidence band
-                    if polygon_points.len() > 2 {
-                        chart.draw_series(std::iter::once(Polygon::new(
-                            polygon_points,
-                            transparent_color.filled(),
-                        )))?;
+                    // Conditionally draw standard deviation bands
+                    if config.show_std_dev {
+                        // Create transparent confidence band
+                        let transparent_color = color.mix(0.3); // 30% opacity
+                        
+                        // Create upper and lower bound points for the filled area
+                        let upper_points: Vec<(f64, f64)> = probability_data
+                            .iter()
+                            .map(|p| (p.evaluations as f64, p.upper_bound))
+                            .collect();
+                        
+                        let mut lower_points: Vec<(f64, f64)> = probability_data
+                            .iter()
+                            .map(|p| (p.evaluations as f64, p.lower_bound))
+                            .collect();
+                        
+                        // Reverse the lower points to create a closed polygon
+                        lower_points.reverse();
+                        
+                        // Combine upper and lower points to form a polygon
+                        let mut polygon_points = upper_points;
+                        polygon_points.extend(lower_points);
+                        
+                        // Draw the filled confidence band
+                        if polygon_points.len() > 2 {
+                            chart.draw_series(std::iter::once(Polygon::new(
+                                polygon_points,
+                                transparent_color.filled(),
+                            )))?;
+                        }
                     }
 
                     // Draw the main line on top of the confidence band
@@ -633,7 +875,11 @@ fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) ->
                         ))?;
 
                     let algo_label = get_algorithm_label(algorithm, config);
-                    let label = format!("{} (avg ± σ)", algo_label);
+                    let label = if config.show_std_dev {
+                        format!("{} (avg ± σ)", algo_label)
+                    } else {
+                        format!("{} (avg)", algo_label)
+                    };
                     series
                         .label(&label)
                         .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
@@ -668,7 +914,11 @@ fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) ->
         }
     }
 
-    chart.configure_series_labels().background_style(&WHITE).border_style(&BLACK).draw()?;
+    chart.configure_series_labels()
+        .background_style(&WHITE)
+        .border_style(&BLACK)
+        .position(SeriesLabelPosition::UpperLeft)
+        .draw()?;
     root.present()?;
 
     println!("Probability plot saved to: {}", output_path.display());