From b0acf3341ffe9fa3b55366eb6dff6a5e1e9a63b7 Mon Sep 17 00:00:00 2001
From: Rutherther <rutherther@ditigal.xyz>
Date: Sun, 2 Nov 2025 20:54:38 +0100
Subject: [PATCH] feat(plotter): add plotting of standard deviation

---
 codes/tsp_plotter/src/main.rs | 212 +++++++++++++++++++++++-----------
 1 file changed, 143 insertions(+), 69 deletions(-)

diff --git a/codes/tsp_plotter/src/main.rs b/codes/tsp_plotter/src/main.rs
index 1054ba4d742c3fcd3a7a0f37257daad8d4ed9c73..0885cc1133c8a91cdcb63aad29644c73e3cd231c 100644
--- a/codes/tsp_plotter/src/main.rs
+++ b/codes/tsp_plotter/src/main.rs
@@ -1,6 +1,7 @@
 use csv::Reader;
 use plotters::prelude::*;
-use std::collections::HashMap;
+use plotters::element::Polygon;
+use std::{collections::HashMap, path::PathBuf};
 use std::fs;
 use std::path::Path;
 use serde::{Deserialize, Serialize};
@@ -26,7 +27,6 @@ struct PlotConfig {
     group_by_algorithm: bool,
     base_path: String,
     output_path: String,
-    optimal_solutions: HashMap<String, f64>,
     targets: Vec<f64>,
     plot_type: PlotType,
     average_targets: bool,
@@ -35,17 +35,12 @@ struct PlotConfig {
 
 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,
@@ -54,6 +49,38 @@ impl Default for PlotConfig {
     }
 }
 
+fn load_optimal_cost(instance_filename: &PathBuf) -> Result<f64, Box<dyn std::error::Error>> {
+    let instance_name = instance_filename
+        .file_stem()
+        .and_then(|s| s.to_str())
+        .ok_or("Could not extract instance name")?
+        .trim_end_matches(".tsp");
+    println!("{:?}", instance_name);
+
+    let solutions_path = instance_filename
+        .parent().unwrap()
+        .parent().unwrap()
+        .parent().unwrap()
+        .join("instances/solutions.txt");
+    println!("{:?}", solutions_path);
+
+    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 calculate_percentage_deviation(fitness: f64, optimal: f64) -> f64 {
     ((fitness - optimal) / optimal) * 100.0
 }
@@ -71,19 +98,19 @@ 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;
@@ -96,7 +123,7 @@ fn create_step_function(data: Vec<DataPoint>) -> Vec<DataPoint> {
             });
         }
     }
-    
+
     result
 }
 
@@ -106,6 +133,15 @@ struct ProbabilityPoint {
     probability: f64,
 }
 
+#[derive(Debug)]
+struct ProbabilityPointWithDeviation {
+    evaluations: u32,
+    probability: f64,
+    std_dev: f64,
+    lower_bound: f64,
+    upper_bound: f64,
+}
+
 fn calculate_success_probability(
     algorithm_data: &HashMap<String, Vec<DataPoint>>,
     target_percentage: f64,
@@ -113,7 +149,7 @@ fn calculate_success_probability(
     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 {
@@ -121,13 +157,13 @@ fn calculate_success_probability(
             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
@@ -136,40 +172,40 @@ fn calculate_success_probability(
                 .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(
+fn calculate_averaged_success_probability_with_deviation(
     algorithm_data: &HashMap<String, Vec<DataPoint>>,
     targets: &[f64],
-) -> Vec<ProbabilityPoint> {
+) -> Vec<ProbabilityPointWithDeviation> {
     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 {
@@ -177,33 +213,46 @@ fn calculate_averaged_success_probability(
             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
+        let mut probabilities = Vec::new();
+
+        // Collect 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;
+                probabilities.push(point.probability);
             }
         }
-        
-        if count > 0 {
-            let averaged_probability = total_probability / count as f64;
-            averaged_points.push(ProbabilityPoint {
+
+        if !probabilities.is_empty() {
+            let mean = probabilities.iter().sum::<f64>() / probabilities.len() as f64;
+            
+            // Calculate standard deviation
+            let variance = probabilities.iter()
+                .map(|p| (p - mean).powi(2))
+                .sum::<f64>() / probabilities.len() as f64;
+            let std_dev = variance.sqrt();
+            
+            // Calculate bounds (mean ± 1 standard deviation, clamped to [0, 1])
+            let lower_bound = (mean - std_dev).max(0.0);
+            let upper_bound = (mean + std_dev).min(1.0);
+            
+            averaged_points.push(ProbabilityPointWithDeviation {
                 evaluations: evaluation,
-                probability: averaged_probability,
+                probability: mean,
+                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)?;
     let mut data = Vec::new();
@@ -222,8 +271,8 @@ fn read_csv_file(file_path: &Path, optimal_solution: f64) -> Result<Vec<DataPoin
     if !data.is_empty() {
         let first_fitness = data[0].fitness;
         let first_percentage = data[0].percentage_deviation;
-        data.insert(0, DataPoint { 
-            evaluations: 1, 
+        data.insert(0, DataPoint {
+            evaluations: 1,
             fitness: first_fitness,
             percentage_deviation: first_percentage
         });
@@ -272,31 +321,26 @@ fn read_plot_data(config: &PlotConfig) -> Result<PlotData, Box<dyn std::error::E
 
     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 optimal_solution = load_optimal_cost(&path).unwrap();
+
                 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);
         }
@@ -307,7 +351,7 @@ fn read_plot_data(config: &PlotConfig) -> Result<PlotData, Box<dyn std::error::E
 
 fn get_color_palette() -> Vec<RGBColor> {
     vec![
-        BLUE, RED, GREEN, CYAN, MAGENTA, 
+        BLUE, RED, GREEN, CYAN, MAGENTA,
         RGBColor(255, 165, 0), // Orange
         RGBColor(128, 0, 128),  // Purple
         RGBColor(255, 192, 203), // Pink
@@ -319,9 +363,9 @@ fn get_color_palette() -> Vec<RGBColor> {
 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)?;
 
@@ -459,7 +503,7 @@ fn create_plot(plot_data: &PlotData, config: &PlotConfig) -> Result<(), Box<dyn
                     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));
@@ -480,9 +524,9 @@ fn create_plot(plot_data: &PlotData, config: &PlotConfig) -> Result<(), Box<dyn
 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)?;
 
@@ -543,15 +587,45 @@ fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) ->
 
     // Plot probability curves
     if config.average_targets {
-        // Plot one averaged curve per algorithm
+        // Plot one averaged curve per algorithm with error bars
         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);
-                
+                let probability_data = calculate_averaged_success_probability_with_deviation(algorithm_data, &config.targets);
+
                 if !probability_data.is_empty() {
                     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(),
+                        )))?;
+                    }
+
+                    // Draw the main line on top of the confidence band
                     let series = chart
                         .draw_series(LineSeries::new(
                             probability_data.iter().map(|p| (p.evaluations as f64, p.probability)),
@@ -559,7 +633,7 @@ 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 = format!("{} (avg ± σ)", algo_label);
                     series
                         .label(&label)
                         .legend(move |(x, y)| PathElement::new(vec![(x, y), (x + 10, y)], &color));
@@ -572,7 +646,7 @@ fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) ->
             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;
@@ -602,7 +676,7 @@ fn create_success_probability_plot(plot_data: &PlotData, config: &PlotConfig) ->
 }
 
 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)?;
@@ -610,28 +684,28 @@ fn load_config(config_path: &str) -> Result<PlotConfig, Box<dyn std::error::Erro
         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);