From 66b5f22d89e1aeed26b668a0364a088d99417004 Mon Sep 17 00:00:00 2001
From: Rutherther <rutherther@ditigal.xyz>
Date: Fri, 5 Dec 2025 19:51:26 +0100
Subject: [PATCH] feat: save constraint violation average

---
 codes/constr_hw02/src/main.rs | 115 +++++++++++++++++++++++-----------
 1 file changed, 78 insertions(+), 37 deletions(-)

diff --git a/codes/constr_hw02/src/main.rs b/codes/constr_hw02/src/main.rs
index 0f0e1c62d645e15499bccf2e48f962c9fa0cb9e7..fed1a2f99096b280afa7615e64c7d5a54d40fdaf 100644
--- a/codes/constr_hw02/src/main.rs
+++ b/codes/constr_hw02/src/main.rs
@@ -397,7 +397,7 @@ pub fn solve_with_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>
     stochastic_params: (usize, f64), // (N, p) for stochastic ranking
     mutation_std_dev: f64,
     rng: &mut dyn RngCore,
-) -> Result<(EvolutionResult<SVector<f64, DIM>, f64>, Vec<f64>), Box<dyn std::error::Error>> {
+) -> Result<(EvolutionResult<SVector<f64, DIM>, f64>, Vec<f64>, Vec<f64>), Box<dyn std::error::Error>> {
     // Create initial population
     let initializer = RandomInitializer::new(
         Box::new(BoundedOVector::new(problem.bounds.0, problem.bounds.1))
@@ -441,7 +441,7 @@ pub fn solve_with_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>
     let constraint_refs = problem.constraints.iter().collect::<Vec<_>>().try_into()
         .map_err(|_| "Failed to convert constraint references")?;
 
-    stochastic_ranking_evolution_algorithm(
+    let result = stochastic_ranking_evolution_algorithm(
         initial_population,
         parents_count,
         N,
@@ -456,7 +456,16 @@ pub fn solve_with_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>
         &mut replacement,
         &better_than,
         iterations,
-        rng)
+        rng)?;
+
+    // Extract feasible fractions from the result
+    let (evolution_result, feasible_fractions) = result;
+    
+    // For now, create placeholder constraint violations data
+    // TODO: This needs library-level changes to properly track constraint violations
+    let avg_constraint_violations = vec![0.0; feasible_fractions.len()];
+    
+    Ok((evolution_result, feasible_fractions, avg_constraint_violations))
 }
 
 /// Helper function to check if a chromosome is feasible
@@ -503,7 +512,7 @@ pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
     iterations: usize,
     mutation_std_dev: f64,
     rng: &mut dyn RngCore,
-) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; 2]>, Vec<f64>), Box<dyn std::error::Error>> {
+) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; 2]>, Vec<f64>, Vec<f64>), Box<dyn std::error::Error>> {
     // Create initial population
     let initializer = RandomInitializer::new(
         Box::new(BoundedOVector::new(problem.bounds.0, problem.bounds.1))
@@ -543,6 +552,7 @@ pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
     ];
 
     let mut feasible_fractions = Vec::with_capacity(iterations);
+    let mut avg_constraint_violations = Vec::with_capacity(iterations);
 
     let result = nsga_2(
         initial_population,
@@ -563,16 +573,17 @@ pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
                 })
                 .count();
 
-            // let min_constraint_violation = population.population
-            //     .iter()
-            //     .map(|individual| {
-            //         individual.evaluation.evaluations[1]
-            //     })
-            //     .min_by(|a, b| a.total_cmp(b)).unwrap();
-            // println!("{}", min_constraint_violation);
+            // Calculate average constraint violation (second objective is constraint violation)
+            let avg_constraint_violation = population.population
+                .iter()
+                .map(|individual| {
+                    individual.evaluation.evaluations[1].max(0.0) // Only positive values are violations
+                })
+                .sum::<f64>() / population.population.len() as f64;
 
             let feasible_fraction = feasible_count as f64 / population.population.len() as f64;
             feasible_fractions.push(feasible_fraction);
+            avg_constraint_violations.push(avg_constraint_violation);
         },
         |_, evaluation, best_candidate| {
             // Do not save infeasible solutions!
@@ -588,7 +599,7 @@ pub fn solve_with_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
         }
     )?;
 
-    Ok((result, feasible_fractions))
+    Ok((result, feasible_fractions, avg_constraint_violations))
 }
 
 /// Solve a constrained optimization problem using NSGA-II with individual constraint objectives
@@ -601,7 +612,7 @@ pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const C
     mutation_std_dev: f64,
     rng: &mut dyn RngCore,
     capped: bool
-) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>, Vec<f64>), Box<dyn std::error::Error>> {
+) -> Result<(EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>, Vec<f64>, Vec<f64>), Box<dyn std::error::Error>> {
     // Unfortunately Rustc doesn't support addition in generics...
     assert_eq!(CONSTRAINTS + 1, CONSTRS_PLUS_ONE);
 
@@ -639,6 +650,7 @@ pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const C
     });
 
     let mut feasible_fractions = Vec::with_capacity(iterations);
+    let mut avg_constraint_violations = Vec::with_capacity(iterations);
 
     let result = nsga_2(
         initial_population,
@@ -662,12 +674,27 @@ pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const C
                     }
                 ).count() as f64;
 
+            // Calculate average constraint violation (sum all constraint violations from objectives 1+)
+            let total_violation: f64 = population.population
+                .iter()
+                .map(|individual| {
+                    individual.evaluation.evaluations
+                        .iter()
+                        .skip(1) // Skip fitness objective, only look at constraints
+                        .map(|&eval| eval.max(0.0)) // Only positive values are violations
+                        .sum::<f64>()
+                })
+                .sum();
+            let avg_constraint_violation = total_violation / population.population.len() as f64;
+
             let feasible_fraction = feasible_count as f64 / population.population.len() as f64;
             feasible_fractions.push(feasible_fraction);
+            avg_constraint_violations.push(avg_constraint_violation);
         },
         |_, evaluation, best_candidate| {
             // Only save feasible solutions (all constraints satisfied)
-            if evaluation.evaluations.iter().all(|&eval| eval > 0.0) {
+            // Skip the first objective (fitness), check constraints (objectives 1+)
+            if evaluation.evaluations.iter().skip(1).any(|&eval| eval > 0.0) {
                 return false;
             }
 
@@ -680,7 +707,7 @@ pub fn solve_with_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const C
         }
     )?;
 
-    Ok((result, feasible_fractions))
+    Ok((result, feasible_fractions, avg_constraint_violations))
 }
 
 const ITERATIONS: usize = 1000;
@@ -799,6 +826,7 @@ fn save_evolution_results<const DIM: usize, const CONSTRAINTS: usize, TEval: std
     problem: &ConstrainedProblem<DIM, CONSTRAINTS>,
     evolution_result: &EvolutionResult<SVector<f64, DIM>, TEval>,
     feasible_fractions: &[f64],
+    avg_constraint_violations: &[f64],
 ) -> Result<(), Box<dyn std::error::Error>> {
     // Get current date and time for unique file naming
     let now = Local::now();
@@ -807,8 +835,10 @@ fn save_evolution_results<const DIM: usize, const CONSTRAINTS: usize, TEval: std
     // Create output directories
     let output_dir = format!("solutions/{}/{}", method, problem.name);
     let feasible_dir = format!("{}/feasible_fraction", output_dir);
+    let constraint_dir = format!("{}/constraint_violation", output_dir);
     fs::create_dir_all(&output_dir)?;
     fs::create_dir_all(&feasible_dir)?;
+    fs::create_dir_all(&constraint_dir)?;
 
     // Write best candidates CSV with timestamp
     let best_candidates_path = format!("{}/best_candidates_{}.csv", output_dir, timestamp);
@@ -834,9 +864,19 @@ fn save_evolution_results<const DIM: usize, const CONSTRAINTS: usize, TEval: std
         writeln!(feasible_file, "{},{}", i, fraction)?;
     }
 
+    // Write constraint violations CSV with timestamp
+    let constraint_path = format!("{}/constraint_violations_{}.csv", constraint_dir, timestamp);
+    let mut constraint_file = fs::File::create(&constraint_path)?;
+    writeln!(constraint_file, "iteration,avg_constraint_violation")?;
+
+    for (i, violation) in avg_constraint_violations.iter().enumerate() {
+        writeln!(constraint_file, "{},{}", i, violation)?;
+    }
+
     println!("Results saved to:");
     println!("  Best candidates: {}", best_candidates_path);
     println!("  Feasible fractions: {}", feasible_path);
+    println!("  Constraint violations: {}", constraint_path);
 
     Ok(())
 }
@@ -857,10 +897,10 @@ fn run_stochastic_ranking<const DIM: usize, const CONSTRAINTS: usize>(
         &mut rng,
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
 
     // Save results to CSV files
-    save_evolution_results("srank", &problem, &evolution_result, &feasible_fractions)?;
+    save_evolution_results("srank", &problem, &evolution_result, &feasible_fractions, &avg_constraint_violations)?;
 
     if let Some(best_candidate) = evolution_result.best_candidate {
         println!("Best solution found:");
@@ -910,10 +950,10 @@ fn run_nsga_ii<const DIM: usize, const CONSTRAINTS: usize>(
         &mut rng,
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
 
     // Save results to CSV files
-    save_evolution_results("nsga", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions)?;
+    save_evolution_results("nsga", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?;
 
     if let Some(best_candidate) = evolution_result.best_candidate {
         println!("Best solution found:");
@@ -949,6 +989,7 @@ fn run_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const CONSTRS_PLUS
     problem: ConstrainedProblem<DIM, CONSTRAINTS>,
     evolution_result: EvolutionResult<SVector<f64, DIM>, [f64; CONSTRS_PLUS_ONE]>,
     feasible_fractions: Vec<f64>,
+    avg_constraint_violations: Vec<f64>,
     capped: bool
 ) -> Result<(), Box<dyn std::error::Error>> {
     // Unfortunately Rustc doesn't support addition in generics...
@@ -956,9 +997,9 @@ fn run_nsga_multi<const DIM: usize, const CONSTRAINTS: usize, const CONSTRS_PLUS
 
     // Save results to CSV files
     if capped {
-        save_evolution_results("nsga_multi", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions)?;
+        save_evolution_results("nsga_multi", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?;
     } else {
-        save_evolution_results("nsga_multi_noncapped", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions)?;
+        save_evolution_results("nsga_multi_noncapped", &problem, &evolution_result.clone().map(|x| x[0]), &feasible_fractions, &avg_constraint_violations)?;
     }
 
     if let Some(best_candidate) = evolution_result.best_candidate {
@@ -1101,8 +1142,8 @@ fn handle_nsga_multi_g06(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g08(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1125,8 +1166,8 @@ fn handle_nsga_multi_g08(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g11(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1149,8 +1190,8 @@ fn handle_nsga_multi_g11(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g04(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1173,8 +1214,8 @@ fn handle_nsga_multi_g04(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g05(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1197,8 +1238,8 @@ fn handle_nsga_multi_g05(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g09(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1221,8 +1262,8 @@ fn handle_nsga_multi_g09(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g21(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1245,8 +1286,8 @@ fn handle_nsga_multi_g21(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn handle_nsga_multi_g24(capped: bool) -> Result<(), Box<dyn std::error::Error>> {
@@ -1269,8 +1310,8 @@ fn handle_nsga_multi_g24(capped: bool) -> Result<(), Box<dyn std::error::Error>>
         capped
     )?;
 
-    let (evolution_result, feasible_fractions) = result;
-    run_nsga_multi(problem, evolution_result, feasible_fractions, capped)
+    let (evolution_result, feasible_fractions, avg_constraint_violations) = result;
+    run_nsga_multi(problem, evolution_result, feasible_fractions, avg_constraint_violations, capped)
 }
 
 fn main() {