changes.diff

diff --git a/benchmarks/pressure_vessel/dependent_space/run.py b/benchmarks/pressure_vessel/dependent_space/run.py
index cf46eaf..b4f9fc7 100644
--- a/benchmarks/pressure_vessel/dependent_space/run.py
+++ b/benchmarks/pressure_vessel/dependent_space/run.py
@@ -139,20 +139,20 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
     π·R²·L + (4/3)·π·R³ ≥ 1,296,000
     → −π·R²·L − (4/3)·π·R³ + 1,296,000 ≤ 0
     """
-    R: float = harmony["R"]
-    Ts: float = harmony["Ts"]
-    Th: float = harmony["Th"]
-    L: float = harmony["L"]
+    r: float = harmony["R"]
+    ts: float = harmony["Ts"]
+    th: float = harmony["Th"]
+    l: float = harmony["L"]
 
     # --- Guard against degenerate designs ---
-    if R <= 0.0 or Ts <= 0.0 or Th <= 0.0 or L <= 0.0:
+    if r <= 0.0 or ts <= 0.0 or th <= 0.0 or l <= 0.0:
         return float("inf"), float("inf")
 
     # --- Fabrication cost ---
-    cost: float = 0.6224 * Ts * R * L + 1.7781 * Th * R**2 + 3.1661 * Ts**2 * L + 19.84 * Ts**2 * R
+    cost: float = 0.6224 * ts * r * l + 1.7781 * th * r**2 + 3.1661 * ts**2 * l + 19.84 * ts**2 * r
 
     # --- Volume constraint (g3) ---
-    volume: float = math.pi * R**2 * L + (4.0 / 3.0) * math.pi * R**3
+    volume: float = math.pi * r**2 * l + (4.0 / 3.0) * math.pi * r**3
     g3: float = -volume + MIN_VOLUME
 
     penalty: float = max(0.0, g3)
@@ -200,7 +200,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Pressure Vessel ΓÇö Dependent Space")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.2f}\n" f"Penalty: {result.best_penalty:.4f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.2f}\nPenalty: {result.best_penalty:.4f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/pressure_vessel/full_parametric_extreme/run.py b/benchmarks/pressure_vessel/full_parametric_extreme/run.py
index aa91858..4fab98f 100644
--- a/benchmarks/pressure_vessel/full_parametric_extreme/run.py
+++ b/benchmarks/pressure_vessel/full_parametric_extreme/run.py
@@ -119,8 +119,8 @@ def _compute_r_min(
         R_min to full machine precision.
     """
 
-    def f(R: float) -> float:
-        return math.pi * R**2 * (l_max + 4.0 * R / 3.0) - volume
+    def f(r: float) -> float:
+        return math.pi * r**2 * (l_max + 4.0 * r / 3.0) - volume
 
     # Bracket: f(0) = -volume < 0;  f(R_upper) >> 0 for any sane problem
     lo, hi = 0.0, R_UPPER
@@ -135,7 +135,7 @@ def _compute_r_min(
     return hi  # conservative (slightly above true root)
 
 
-def _l_min_from_volume(R: float, volume: float) -> float:
+def _l_min_from_volume(r: float, volume: float) -> float:
     """
     Invert the volume equation for the minimum cylindrical length.
 
@@ -156,11 +156,11 @@ def _l_min_from_volume(R: float, volume: float) -> float:
     float
         Lower bound for L.
     """
-    head_volume: float = (4.0 / 3.0) * math.pi * R**3
+    head_volume: float = (4.0 / 3.0) * math.pi * r**3
     remaining: float = volume - head_volume
     if remaining <= 0.0:
         return L_LOWER  # heads alone satisfy volume
-    return max(L_LOWER, remaining / (math.pi * R**2))
+    return max(L_LOWER, remaining / (math.pi * r**2))
 
 
 # ---------------------------------------------------------------------------
@@ -233,12 +233,17 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
     3.1661 ┬╖ Ts┬▓ ┬╖ L       ΓÇö longitudinal shell forming
     19.84  ┬╖ Ts┬▓ ┬╖ R       ΓÇö circumferential shell forming
     """
-    R: float = harmony["R"]
-    Ts: float = harmony["Ts"]
-    Th: float = harmony["Th"]
-    L: float = harmony["L"]
-
-    cost: float = 0.6224 * Ts * R * L + 1.7781 * Th * R**2 + 3.1661 * Ts**2 * L + 19.84 * Ts**2 * R
+    r: float = harmony["R"]
+    ts: float = harmony["Ts"]
+    th: float = harmony["Th"]
+    l: float = harmony["L"]
+
+    cost: float = (
+        0.6224 * ts * r * l
+        + 1.7781 * th * r**2
+        + 3.1661 * ts**2 * l
+        + 19.84 * ts**2 * r
+    )
 
     return cost, 0.0
 
@@ -287,10 +292,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Pressure Vessel ΓÇö Full Parametric Extreme")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.2f}\n"
-        f"Penalty: {result.best_penalty:.4f}\n"
-        f"R_min: {R_MIN:.2f} in\n"
-        f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.2f}\nPenalty: {result.best_penalty:.4f}\nR_min: {R_MIN:.2f} in\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/pressure_vessel/semi_dependent/run.py b/benchmarks/pressure_vessel/semi_dependent/run.py
index 17d538d..4cc8252 100644
--- a/benchmarks/pressure_vessel/semi_dependent/run.py
+++ b/benchmarks/pressure_vessel/semi_dependent/run.py
@@ -73,7 +73,7 @@ HEAVY_PENALTY: float = 1e9  # Penalty for geometrically infeasible R
 # ---------------------------------------------------------------------------
 # Helper: volume-derived minimum L
 # ---------------------------------------------------------------------------
-def _l_min_from_volume(R: float) -> float:
+def _l_min_from_volume(r: float) -> float:
     """
     Invert the volume equation to find the minimum cylindrical length
     that satisfies g3 for a given inner radius R.
@@ -85,11 +85,11 @@ def _l_min_from_volume(R: float) -> float:
     sufficient volume, making L_min negative.  The result is clamped
     to a minimum of 10.0 (the benchmark lower bound for L).
     """
-    head_volume: float = (4.0 / 3.0) * math.pi * R**3
+    head_volume: float = (4.0 / 3.0) * math.pi * r**3
     remaining: float = MIN_VOLUME - head_volume
     if remaining <= 0:
         return 10.0  # heads alone satisfy volume
-    return remaining / (math.pi * R**2)
+    return remaining / (math.pi * r**2)
 
 
 # ---------------------------------------------------------------------------
@@ -160,22 +160,22 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
     For all other designs the penalty is 0.0 because the space
     guarantees g1ΓÇôg4 satisfaction.
     """
-    R: float = harmony["R"]
-    Ts: float = harmony["Ts"]
-    Th: float = harmony["Th"]
-    L: float = harmony["L"]
+    r: float = harmony["R"]
+    ts: float = harmony["Ts"]
+    th: float = harmony["Th"]
+    l: float = harmony["L"]
 
-    if R <= 0.0 or Ts <= 0.0 or Th <= 0.0 or L <= 0.0:
+    if r <= 0.0 or ts <= 0.0 or th <= 0.0 or l <= 0.0:
         return float("inf"), float("inf")
 
     # --- Check geometric feasibility of R ---------------------------------
-    l_min_required: float = _l_min_from_volume(R)
+    l_min_required: float = _l_min_from_volume(r)
     if l_min_required > 240.0:
         # This R is too small ΓÇö no valid L can satisfy the volume.
         return HEAVY_PENALTY, HEAVY_PENALTY
 
     # --- Fabrication cost --------------------------------------------------
-    cost: float = 0.6224 * Ts * R * L + 1.7781 * Th * R**2 + 3.1661 * Ts**2 * L + 19.84 * Ts**2 * R
+    cost: float = 0.6224 * ts * r * l + 1.7781 * th * r**2 + 3.1661 * ts**2 * l + 19.84 * ts**2 * r
 
     return cost, 0.0
 
@@ -220,7 +220,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Pressure Vessel ΓÇö Semi-Dependent Space")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.2f}\n" f"Penalty: {result.best_penalty:.4f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.2f}\nPenalty: {result.best_penalty:.4f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/pressure_vessel/static_penalty/run.py b/benchmarks/pressure_vessel/static_penalty/run.py
index 4b9174d..4ff1990 100644
--- a/benchmarks/pressure_vessel/static_penalty/run.py
+++ b/benchmarks/pressure_vessel/static_penalty/run.py
@@ -82,23 +82,23 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
     Term 3:  3.1661 ┬╖ Ts┬▓ ┬╖ L         ΓÇö longitudinal shell forming
     Term 4: 19.84   ┬╖ Ts┬▓ ┬╖ R         ΓÇö circumferential shell forming
     """
-    Ts: float = harmony["Ts"]
-    Th: float = harmony["Th"]
-    R: float = harmony["R"]
-    L: float = harmony["L"]
+    ts: float = harmony["Ts"]
+    th: float = harmony["Th"]
+    r: float = harmony["R"]
+    l: float = harmony["L"]
 
     # --- Guard against degenerate designs ---
-    if Ts <= 0 or Th <= 0 or R <= 0 or L <= 0:
+    if ts <= 0 or th <= 0 or r <= 0 or l <= 0:
         return float("inf"), float("inf")
 
     # --- Fabrication cost ---
-    cost: float = 0.6224 * Ts * R * L + 1.7781 * Th * R**2 + 3.1661 * Ts**2 * L + 19.84 * Ts**2 * R
+    cost: float = 0.6224 * ts * r * l + 1.7781 * th * r**2 + 3.1661 * ts**2 * l + 19.84 * ts**2 * r
 
     # --- All 4 constraints ---
-    g1 = -Ts + 0.0193 * R  # shell hoop
-    g2 = -Th + 0.00954 * R  # head hoop
-    g3 = -(math.pi * R**2 * L) - (4.0 / 3.0) * math.pi * R**3 + MIN_VOLUME  # volume
-    g4 = L - 240.0  # length
+    g1 = -ts + 0.0193 * r  # shell hoop
+    g2 = -th + 0.00954 * r  # head hoop
+    g3 = -(math.pi * r**2 * l) - (4.0 / 3.0) * math.pi * r**3 + MIN_VOLUME  # volume
+    g4 = l - 240.0  # length
 
     penalty: float = sum(max(0.0, g) for g in [g1, g2, g3, g4])
 
@@ -145,7 +145,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Pressure Vessel ΓÇö Static Penalty")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.2f}\n" f"Penalty: {result.best_penalty:.4f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.2f}\nPenalty: {result.best_penalty:.4f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/retaining_wall/dependent_space/run.py b/benchmarks/retaining_wall/dependent_space/run.py
index c6a2613..50bf2ed 100644
--- a/benchmarks/retaining_wall/dependent_space/run.py
+++ b/benchmarks/retaining_wall/dependent_space/run.py
@@ -60,10 +60,10 @@ def _x3_min(ctx: Dict[str, float]) -> float:
     fc = round(ctx["fc"] / 5.0) * 5.0
     x4 = round(ctx["x4"] / 50.0) * 50.0
 
-    Vu_stem = 1.6 * (0.5 * KA * GAMMA_SOIL * H_STEM**2) * 1000.0  # N
+    vu_stem = 1.6 * (0.5 * KA * GAMMA_SOIL * H_STEM**2) * 1000.0  # N
     # phi * Vc = 0.75 * 0.17 * sqrt(f'c) * bw * d
     shear_cap_per_mm = 0.75 * 0.17 * math.sqrt(fc) * WIDTH_B
-    d_req = Vu_stem / shear_cap_per_mm
+    d_req = vu_stem / shear_cap_per_mm
 
     return float(max(300.0, x4, math.ceil((d_req + COVER) / 50.0) * 50.0))
 
@@ -75,8 +75,8 @@ def check_geotech_stability(x1, x2, x3, x4, x5) -> bool:
     if l_heel < 0:
         return False
 
-    Pa = 0.5 * KA * GAMMA_SOIL * (h_tot**2)
-    Mo = Pa * (h_tot / 3.0)
+    pa = 0.5 * KA * GAMMA_SOIL * (h_tot**2)
+    mo = pa * (h_tot / 3.0)
 
     w_base = (x1 / 1000.0) * t_base_m * GAMMA_CONC
     w_base_x = (x1 / 1000.0) / 2.0
@@ -90,12 +90,12 @@ def check_geotech_stability(x1, x2, x3, x4, x5) -> bool:
     w_soil = (l_heel / 1000.0) * H_STEM * GAMMA_SOIL
     w_soil_x = (x1 - l_heel / 2.0) / 1000.0
 
-    sum_W = w_base + w_stem_rect + w_stem_tri + w_soil
-    Mr = (w_base * w_base_x) + (w_stem_rect * w_stem_rect_x) + (w_stem_tri * w_stem_tri_x) + (w_soil * w_soil_x)
+    sum_w = w_base + w_stem_rect + w_stem_tri + w_soil
+    mr = (w_base * w_base_x) + (w_stem_rect * w_stem_rect_x) + (w_stem_tri * w_stem_tri_x) + (w_soil * w_soil_x)
 
-    fs_ov = Mr / Mo if Mo > 0 else 100.0
-    fs_sl = (sum_W * math.tan(PHI_SOIL)) / Pa
-    e = abs((x1 / 2000.0) - ((Mr - Mo) / sum_W))
+    fs_ov = mr / mo if mo > 0 else 100.0
+    fs_sl = (sum_w * math.tan(PHI_SOIL)) / pa
+    e = abs((x1 / 2000.0) - ((mr - mo) / sum_w))
 
     return fs_ov >= 2.5 and fs_sl >= 2.5 and e <= (x1 / 6000.0)
 
@@ -171,18 +171,18 @@ def objective(config: Dict[str, Any]) -> Tuple[float, float]:
     l_heel = x1 - x2 - x3
 
     # Demands
-    Pa_stem = 0.5 * KA * GAMMA_SOIL * (H_STEM**2)
-    Mu_stem = 1.6 * Pa_stem * (H_STEM / 3.0)
+    pa_stem = 0.5 * KA * GAMMA_SOIL * (H_STEM**2)
+    mu_stem = 1.6 * pa_stem * (H_STEM / 3.0)
     d_stem = x3 - COVER
 
     w_heel_total = 1.2 * (GAMMA_SOIL * H_STEM + GAMMA_CONC * t_base_m)
-    Vu_heel_N = w_heel_total * (l_heel / 1000.0) * 1000.0
-    Mu_heel = w_heel_total * ((l_heel / 1000.0) ** 2 / 2.0)
+    vu_heel_n = w_heel_total * (l_heel / 1000.0) * 1000.0
+    mu_heel = w_heel_total * ((l_heel / 1000.0) ** 2 / 2.0)
     d_base = x5 - COVER
 
     # Embed shear check for heel (dynamically penalized if it fails because x5 wasn't pruned)
-    phi_Vc_heel = 0.75 * 0.17 * math.sqrt(fc) * WIDTH_B * d_base
-    shear_pen = max(0.0, Vu_heel_N - phi_Vc_heel) ** 2 * 1e4
+    phi_vc_heel = 0.75 * 0.17 * math.sqrt(fc) * WIDTH_B * d_base
+    shear_pen = max(0.0, vu_heel_n - phi_vc_heel) ** 2 * 1e4
 
     code_stem = config["dc_stem"]
     code_base = config["dc_base"]
@@ -190,21 +190,21 @@ def objective(config: Dict[str, Any]) -> Tuple[float, float]:
     if code_stem is None or code_base is None or l_heel < 0:
         return 1e12, 1e8 + space_penalty
 
-    As_s = steel_area(code_stem)
-    As_b = steel_area(code_base)
+    as_s = steel_area(code_stem)
+    as_b = steel_area(code_base)
 
     # Check Flexure (Because ACIRebar handles rho & fit, but we must verify Moment limit)
-    a_stem = (As_s * FY_MPA) / (0.85 * fc * WIDTH_B)
-    phi_Mn_stem = 0.9 * As_s * FY_MPA * max(0.1, d_stem - a_stem / 2.0) * 1e-6
-    mu_pen_stem = max(0.0, Mu_stem - phi_Mn_stem) ** 2 * 1e4
+    a_stem = (as_s * FY_MPA) / (0.85 * fc * WIDTH_B)
+    phi_mn_stem = 0.9 * as_s * FY_MPA * max(0.1, d_stem - a_stem / 2.0) * 1e-6
+    mu_pen_stem = max(0.0, mu_stem - phi_mn_stem) ** 2 * 1e4
 
-    a_base = (As_b * FY_MPA) / (0.85 * fc * WIDTH_B)
-    phi_Mn_base = 0.9 * As_b * FY_MPA * max(0.1, d_base - a_base / 2.0) * 1e-6
-    mu_pen_base = max(0.0, Mu_heel - phi_Mn_base) ** 2 * 1e4
+    a_base = (as_b * FY_MPA) / (0.85 * fc * WIDTH_B)
+    phi_mn_base = 0.9 * as_b * FY_MPA * max(0.1, d_base - a_base / 2.0) * 1e-6
+    mu_pen_base = max(0.0, mu_heel - phi_mn_base) ** 2 * 1e4
 
     vol_conc = (x1 / 1000.0) * t_base_m + 0.5 * (x3 + x4) / 1000.0 * H_STEM
     cost_conc = (0.5 + 0.02 * fc) * vol_conc
-    cost_steel = 50.0 * (As_s + As_b) / 1000.0
+    cost_steel = 50.0 * (as_s + as_b) / 1000.0
     total_cost = cost_conc + cost_steel
 
     return total_cost, space_penalty + shear_pen + mu_pen_stem + mu_pen_base
diff --git a/benchmarks/retaining_wall/static_penalty/run.py b/benchmarks/retaining_wall/static_penalty/run.py
index 084f4c6..9452418 100644
--- a/benchmarks/retaining_wall/static_penalty/run.py
+++ b/benchmarks/retaining_wall/static_penalty/run.py
@@ -83,8 +83,8 @@ def objective(config: Dict[str, Any]) -> Tuple[float, float]:
     h_tot = H_STEM + t_base_m
 
     # Active Earth Thrust
-    Pa = 0.5 * KA * GAMMA_SOIL * (h_tot**2)
-    Mo = Pa * (h_tot / 3.0)
+    pa = 0.5 * KA * GAMMA_SOIL * (h_tot**2)
+    mo = pa * (h_tot / 3.0)
 
     # Weights and Centers
     w_base = (x1 / 1000.0) * t_base_m * GAMMA_CONC
@@ -99,61 +99,61 @@ def objective(config: Dict[str, Any]) -> Tuple[float, float]:
     w_soil = (l_heel / 1000.0) * H_STEM * GAMMA_SOIL
     w_soil_x = (x1 - l_heel / 2.0) / 1000.0
 
-    sum_W = w_base + w_stem_rect + w_stem_tri + w_soil
-    Mr = (w_base * w_base_x) + (w_stem_rect * w_stem_rect_x) + (w_stem_tri * w_stem_tri_x) + (w_soil * w_soil_x)
+    sum_w = w_base + w_stem_rect + w_stem_tri + w_soil
+    mr = (w_base * w_base_x) + (w_stem_rect * w_stem_rect_x) + (w_stem_tri * w_stem_tri_x) + (w_soil * w_soil_x)
 
-    fs_overturning = Mr / Mo if Mo > 0 else 100.0
-    fs_sliding = (sum_W * math.tan(PHI_SOIL)) / Pa
+    fs_overturning = mr / mo if mo > 0 else 100.0
+    fs_sliding = (sum_w * math.tan(PHI_SOIL)) / pa
 
     # Eccentricity
-    e = (x1 / 2000.0) - ((Mr - Mo) / sum_W)
+    e = (x1 / 2000.0) - ((mr - mo) / sum_w)
 
     # --- Structural ACI Demands ---
     # Stem at base
-    Pa_stem = 0.5 * KA * GAMMA_SOIL * (H_STEM**2)
-    Vu_stem = 1.6 * Pa_stem
-    Mu_stem = 1.6 * Pa_stem * (H_STEM / 3.0)
+    pa_stem = 0.5 * KA * GAMMA_SOIL * (H_STEM**2)
+    vu_stem = 1.6 * pa_stem
+    mu_stem = 1.6 * pa_stem * (H_STEM / 3.0)
     d_stem = x3 - COVER
 
     # Heel at back of stem
     w_heel_total = 1.2 * (GAMMA_SOIL * H_STEM + GAMMA_CONC * t_base_m)
-    Vu_heel = w_heel_total * (l_heel / 1000.0)
-    Mu_heel = w_heel_total * ((l_heel / 1000.0) ** 2 / 2.0)
+    vu_heel = w_heel_total * (l_heel / 1000.0)
+    mu_heel = w_heel_total * ((l_heel / 1000.0) ** 2 / 2.0)
     d_base = x5 - COVER
 
     # Rebars
     db_s = DIAMS[idx_stem_dia]
-    As_s = float(n_s) * math.pi * (db_s**2) / 4.0
+    as_s = float(n_s) * math.pi * (db_s**2) / 4.0
 
     db_b = DIAMS[idx_base_dia]
-    As_b = float(n_b) * math.pi * (db_b**2) / 4.0
+    as_b = float(n_b) * math.pi * (db_b**2) / 4.0
 
     # --- Structural Penalties (g_i <= 0) ---
     beta1 = get_beta1(fc)
     rho_min = max(0.25 * math.sqrt(fc) / FY_MPA, 1.4 / FY_MPA)
     rho_max = 0.85 * beta1 * (fc / FY_MPA) * (3.0 / 8.0)
 
-    def check_aci(M_u, V_u, d_eff, A_s, d_b, n):
-        phi_Vc = 0.75 * 0.17 * math.sqrt(fc) * WIDTH_B * (d_eff / 1000.0)
-        a = (A_s * FY_MPA) / (0.85 * fc * WIDTH_B)
-        phi_Mn = 0.9 * A_s * FY_MPA * max(0.1, d_eff - a / 2.0) * 1e-6
-        rho = A_s / (WIDTH_B * d_eff)
+    def check_aci(m_u, v_u, d_eff, a_s, d_b, n):
+        phi_vc = 0.75 * 0.17 * math.sqrt(fc) * WIDTH_B * (d_eff / 1000.0)
+        a = (a_s * FY_MPA) / (0.85 * fc * WIDTH_B)
+        phi_mn = 0.9 * a_s * FY_MPA * max(0.1, d_eff - a / 2.0) * 1e-6
+        rho = a_s / (WIDTH_B * d_eff)
         # Spacing
         s_available = (WIDTH_B - 2.0 * COVER - n * d_b) / max(1, n - 1)
         s_min = max(25.0, d_b)
 
-        return [M_u - phi_Mn, V_u - phi_Vc, rho_min - rho, rho - rho_max, s_min - s_available]
+        return [m_u - phi_mn, v_u - phi_vc, rho_min - rho, rho - rho_max, s_min - s_available]
 
     g = [2.5 - fs_overturning, 2.5 - fs_sliding, abs(e) - (x1 / 6000.0)]
-    g.extend(check_aci(Mu_stem, Vu_stem, d_stem, As_s, db_s, n_s))
-    g.extend(check_aci(Mu_heel, Vu_heel, d_base, As_b, db_b, n_b))
+    g.extend(check_aci(mu_stem, vu_stem, d_stem, as_s, db_s, n_s))
+    g.extend(check_aci(mu_heel, vu_heel, d_base, as_b, db_b, n_b))
 
     penalty_value = sum(max(0.0, limit) ** 2 for limit in g) * 1e6
 
     # --- Cost Evaluation ---
     vol_conc = (x1 / 1000.0) * t_base_m + 0.5 * (x3 + x4) / 1000.0 * H_STEM
     cost_conc = (0.5 + 0.02 * fc) * vol_conc
-    cost_steel = 50.0 * (As_s + As_b) / 1000.0
+    cost_steel = 50.0 * (as_s + as_b) / 1000.0
     total_cost = cost_conc + cost_steel
 
     return total_cost, penalty_value
diff --git a/benchmarks/spring_design/semi_dependent/run.py b/benchmarks/spring_design/semi_dependent/run.py
index d17918c..db61ef4 100644
--- a/benchmarks/spring_design/semi_dependent/run.py
+++ b/benchmarks/spring_design/semi_dependent/run.py
@@ -183,7 +183,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Spring Design ΓÇö Semi-Dependent Space")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.6f}\n" f"Penalty: {result.best_penalty:.6f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.6f}\nPenalty: {result.best_penalty:.6f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/spring_design/static_penalty/run.py b/benchmarks/spring_design/static_penalty/run.py
index 6105d65..1a9a97c 100644
--- a/benchmarks/spring_design/static_penalty/run.py
+++ b/benchmarks/spring_design/static_penalty/run.py
@@ -141,7 +141,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Spring Design ΓÇö Static Penalty")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.6f}\n" f"Penalty: {result.best_penalty:.6f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.6f}\nPenalty: {result.best_penalty:.6f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/utils/plotter.py b/benchmarks/utils/plotter.py
index cf3b658..59d2bae 100644
--- a/benchmarks/utils/plotter.py
+++ b/benchmarks/utils/plotter.py
@@ -225,30 +225,8 @@ class ConvergencePlotter:
         return self
 
     # --------------------------------------------------------------- render
-    def plot(self, save_path: Optional[Union[str, Path]] = None) -> None:
-        """
-        Render the convergence plot.
-
-        Three scenarios are handled automatically:
-
-        1. **All penalties zero** ΓÇö only the fitness curve is drawn.
-        2. **No feasible solution found** ΓÇö only the penalty curve is drawn.
-        3. **Mixed** ΓÇö penalty is drawn (shifted) up to the feasibility
-           transition point, then fitness takes over.
-
-        Parameters
-        ----------
-        save_path : str | Path | None
-            If given, the figure is saved at 300 DPI.  Otherwise ``plt.show()``
-            is called.
-        """
-        _apply_academic_style()
-
-        if not self.iterations_data:
-            print("No data to plot!")
-            return
-
-        # --- Phase detection ------------------------------------------------
+    def _phase_detection(self) -> Tuple[int, float, bool]:
+        """Detect the first iteration where penalty reaches (close to) zero."""
         first_zero_idx: int = -1
         first_feasible_fitness: float = -1.0
 
@@ -259,128 +237,195 @@ class ConvergencePlotter:
                 break
 
         all_zero: bool = all(np.isclose(r.penalty, 0.0) for r in self.iterations_data)
+        return first_zero_idx, first_feasible_fitness, all_zero
 
+    def _extract_series(
+        self,
+    ) -> Tuple[
+        List[int],
+        List[float],
+        List[float],
+        Optional[float],
+        Optional[float],
+        float,
+        Optional[float],
+    ]:
         iterations = [r.iteration for r in self.iterations_data]
         fitness_values = [r.fitness for r in self.iterations_data]
         penalty_values = [r.penalty for r in self.iterations_data]
 
-        # Pre-compute stats
         min_fitness: Optional[float] = min(fitness_values) if fitness_values else None
         last_fitness: Optional[float] = fitness_values[-1] if fitness_values else None
         max_penalty: float = max(penalty_values) if penalty_values else 0.0
         first_fitness: Optional[float] = fitness_values[0] if fitness_values else None
 
-        # Physical y-values for reference lines
-        y_pen_zero_phys: Optional[float] = 0.0
-        y_pen_max_phys: Optional[float] = max_penalty
+        return (
+            iterations,
+            fitness_values,
+            penalty_values,
+            min_fitness,
+            last_fitness,
+            max_penalty,
+            first_fitness,
+        )
+
+    def _build_draw_series(
+        self,
+        *,
+        all_zero: bool,
+        first_zero_idx: int,
+        first_feasible_fitness: float,
+        fitness_values: List[float],
+        penalty_values: List[float],
+        min_fitness: Optional[float],
+        last_fitness: Optional[float],
+        max_penalty: float,
+    ) -> Tuple[List[Optional[float]], List[Optional[float]], Optional[float], Optional[float], Optional[float], Optional[float]]:
+        """
+        Build the two draw-series arrays (with None gaps) and update axis/reference stats.
 
-        # --- Build draw series (penaltyΓÇôfitness alignment) ------------------
+        Returns:
+            penalties_draw, fitness_draw, min_fitness, last_fitness, y_pen_zero_phys, y_pen_max_phys
+        """
         penalties_draw: List[Optional[float]] = []
         fitness_draw: List[Optional[float]] = []
 
+        y_pen_zero_phys: Optional[float] = 0.0
+        y_pen_max_phys: Optional[float] = max_penalty
+
         if all_zero:
-            # Scenario 1: purely feasible run ΓÇö no penalty curve
             penalties_draw = [None] * len(self.iterations_data)
             fitness_draw = list(fitness_values)
             y_pen_zero_phys = None
             y_pen_max_phys = None
+            return penalties_draw, fitness_draw, min_fitness, last_fitness, y_pen_zero_phys, y_pen_max_phys
 
-        elif first_zero_idx == -1:
-            # Scenario 2: never reached feasibility ΓÇö no fitness curve
+        if first_zero_idx == -1:
             penalties_draw = list(penalty_values)
             fitness_draw = [None] * len(self.iterations_data)
-            min_fitness = None
-            last_fitness = None
             y_pen_zero_phys = None
+            y_pen_max_phys = None
+            return penalties_draw, fitness_draw, None, None, y_pen_zero_phys, y_pen_max_phys
 
-        else:
-            # Scenario 3: transition ΓÇö shift penalties up by first feasible
-            for i, rec in enumerate(self.iterations_data):
-                if i <= first_zero_idx:
-                    penalties_draw.append(rec.penalty + first_feasible_fitness)
-                    fitness_draw.append(None)
-                else:
-                    penalties_draw.append(None)
-                    fitness_draw.append(rec.fitness)
-
-            y_pen_zero_phys = 0.0 + first_feasible_fitness
-            y_pen_max_phys = max_penalty + first_feasible_fitness
-            min_fitness = last_fitness
-
-        # --- Figure ---------------------------------------------------------
-        fig, ax1 = plt.subplots(figsize=(12, 6))
-        ax1.set_xlabel(self.x_label)
+        for i, rec in enumerate(self.iterations_data):
+            if i <= first_zero_idx:
+                penalties_draw.append(rec.penalty + first_feasible_fitness)
+                fitness_draw.append(None)
+            else:
+                penalties_draw.append(None)
+                fitness_draw.append(rec.fitness)
 
-        valid_fit_iters = [iterations[i] for i, v in enumerate(fitness_draw) if v is not None]
-        valid_fit_vals = [v for v in fitness_draw if v is not None]
+        y_pen_zero_phys = 0.0 + first_feasible_fitness
+        y_pen_max_phys = max_penalty + first_feasible_fitness
+        min_fitness = last_fitness  # match original plotting logic
+        return penalties_draw, fitness_draw, min_fitness, last_fitness, y_pen_zero_phys, y_pen_max_phys
 
-        if valid_fit_vals:
-            ax1.set_ylabel(self.y_fit_label, color="#2ca02c")
-            ax1.tick_params(axis="y", labelcolor="#2ca02c")
-            ax1.plot(
-                valid_fit_iters,
-                valid_fit_vals,
-                color="#2ca02c",
-                linestyle="-",
-                label=self.y_fit_label,
-                linewidth=2,
-            )
-        else:
+    def _plot_fitness_curve(
+        self,
+        *,
+        ax1: plt.Axes,
+        iterations: List[int],
+        fitness_draw: List[Optional[float]],
+    ) -> List[float]:
+        valid_fit_vals = [v for v in fitness_draw if v is not None]
+        if not valid_fit_vals:
             ax1.set_ylabel("")
             ax1.set_yticks([])
+            return []
 
-        ax2 = ax1.twinx()
-        if not valid_fit_vals:
-            ax2.grid(True, linestyle="--", alpha=0.55)
+        valid_fit_iters = [iterations[i] for i, v in enumerate(fitness_draw) if v is not None]
+        ax1.set_ylabel(self.y_fit_label, color="#2ca02c")
+        ax1.tick_params(axis="y", labelcolor="#2ca02c")
+        ax1.plot(
+            valid_fit_iters,
+            valid_fit_vals,
+            color="#2ca02c",
+            linestyle="-",
+            label=self.y_fit_label,
+            linewidth=2,
+        )
+        return valid_fit_vals
 
-        valid_pen_iters = [iterations[i] for i, v in enumerate(penalties_draw) if v is not None]
+    def _plot_penalty_curve(
+        self,
+        *,
+        ax2: plt.Axes,
+        iterations: List[int],
+        penalties_draw: List[Optional[float]],
+        valid_fit_vals: List[float],
+    ) -> List[float]:
         valid_pen_vals = [v for v in penalties_draw if v is not None]
+        if not valid_fit_vals:
+            ax2.grid(True, linestyle="--", alpha=0.55)
 
-        if valid_pen_vals:
-            ax2.set_ylabel(self.y_pen_label, color="#d62728")
-            ax2.tick_params(axis="y", labelcolor="#d62728")
-            ax2.plot(
-                valid_pen_iters,
-                valid_pen_vals,
-                color="#d62728",
-                linestyle="-",
-                label=self.y_pen_label,
-                linewidth=2,
-            )
-        else:
+        if not valid_pen_vals:
             ax2.set_ylabel("")
             ax2.set_yticks([])
+            return []
 
-        # --- Axis scaling ---------------------------------------------------
+        valid_pen_iters = [iterations[i] for i, v in enumerate(penalties_draw) if v is not None]
+        ax2.set_ylabel(self.y_pen_label, color="#d62728")
+        ax2.tick_params(axis="y", labelcolor="#d62728")
+        ax2.plot(
+            valid_pen_iters,
+            valid_pen_vals,
+            color="#d62728",
+            linestyle="-",
+            label=self.y_pen_label,
+            linewidth=2,
+        )
+        return valid_pen_vals
+
+    def _apply_axis_scaling(
+        self,
+        *,
+        ax1: plt.Axes,
+        ax2: plt.Axes,
+        valid_fit_vals: List[float],
+        valid_pen_vals: List[float],
+        min_fitness: Optional[float],
+        first_fitness: Optional[float],
+        all_zero: bool,
+        y_pen_zero_phys: Optional[float],
+        y_pen_max_phys: Optional[float],
+    ) -> None:
         if self.y_lim is not None:
             ax1.set_ylim(*self.y_lim)
             ax2.set_ylim(*self.y_lim)
-        else:
-            all_y: List[float] = []
-            if valid_fit_vals:
-                all_y.extend(valid_fit_vals)
-            if valid_pen_vals:
-                all_y.extend(valid_pen_vals)
-            if min_fitness is not None:
-                all_y.append(min_fitness)
-            if all_zero and first_fitness is not None:
-                all_y.append(first_fitness)
-            if y_pen_zero_phys is not None:
-                all_y.append(y_pen_zero_phys)
-            if y_pen_max_phys is not None:
-                all_y.append(y_pen_max_phys)
-
-            if all_y:
-                lo, hi = min(all_y), max(all_y)
-                pad = 0.5 if hi == lo else (hi - lo) * 0.05
-                ax1.set_ylim(lo - pad, hi + pad)
-                ax2.set_ylim(lo - pad, hi + pad)
+            return
 
-        if self.x_lim is not None:
-            ax1.set_xlim(*self.x_lim)
-            ax2.set_xlim(*self.x_lim)
+        all_y: List[float] = []
+        if valid_fit_vals:
+            all_y.extend(valid_fit_vals)
+        if valid_pen_vals:
+            all_y.extend(valid_pen_vals)
+        if min_fitness is not None:
+            all_y.append(min_fitness)
+        if all_zero and first_fitness is not None:
+            all_y.append(first_fitness)
+        if y_pen_zero_phys is not None:
+            all_y.append(y_pen_zero_phys)
+        if y_pen_max_phys is not None:
+            all_y.append(y_pen_max_phys)
+
+        if not all_y:
+            return
+
+        lo, hi = min(all_y), max(all_y)
+        pad = 0.5 if hi == lo else (hi - lo) * 0.05
+        ax1.set_ylim(lo - pad, hi + pad)
+        ax2.set_ylim(lo - pad, hi + pad)
 
-        # --- Tick formatters (penaltyΓÇôfitness alignment) --------------------
+    def _configure_tick_formatters(
+        self,
+        *,
+        ax1: plt.Axes,
+        ax2: plt.Axes,
+        first_zero_idx: int,
+        first_feasible_fitness: float,
+        all_zero: bool,
+        valid_fit_vals: List[float],
+    ) -> None:
         _fzidx = first_zero_idx
         _fffit = first_feasible_fitness if first_zero_idx != -1 else None
         _allz = all_zero
@@ -403,7 +448,18 @@ class ConvergencePlotter:
         if valid_fit_vals:
             ax1.yaxis.set_major_formatter(mticker.FuncFormatter(_fitness_formatter))
 
-        # --- Reference lines ------------------------------------------------
+    def _draw_reference_lines(
+        self,
+        *,
+        ax1: plt.Axes,
+        ax2: plt.Axes,
+        first_zero_idx: int,
+        all_zero: bool,
+        first_feasible_fitness: float,
+        min_fitness: Optional[float],
+        last_fitness: Optional[float],
+        y_pen_zero_phys: Optional[float],
+    ) -> None:
         if first_zero_idx != -1 and not all_zero and first_feasible_fitness is not None:
             ax1.axhline(
                 first_feasible_fitness,
@@ -444,23 +500,22 @@ class ConvergencePlotter:
                     alpha=0.6,
                 )
 
-        ax1.set_title(self.title, pad=15)
-
-        # --- Unified legend -------------------------------------------------
+    def _draw_legend(self, *, ax1: plt.Axes, ax2: plt.Axes) -> None:
         lines1, labels1 = ax1.get_legend_handles_labels()
         lines2, labels2 = ax2.get_legend_handles_labels()
         all_lines = lines1 + lines2
+        if not all_lines:
+            return
         all_labels = labels1 + labels2
-        if all_lines:
-            ax2.legend(
-                all_lines,
-                all_labels,
-                loc="upper right",
-                framealpha=0.85,
-                edgecolor="#cccccc",
-            )
+        ax2.legend(
+            all_lines,
+            all_labels,
+            loc="upper right",
+            framealpha=0.85,
+            edgecolor="#cccccc",
+        )
 
-        # --- User-defined info boxes ----------------------------------------
+    def _draw_info_boxes(self, *, fig: plt.Figure, ax2: plt.Axes) -> None:
         for box in self.text_boxes:
             fig.text(
                 box.x,
@@ -479,8 +534,113 @@ class ConvergencePlotter:
                 ),
             )
 
-        fig.tight_layout()
+    def plot(self, save_path: Optional[Union[str, Path]] = None) -> None:
+        """
+        Render the convergence plot.
+
+        Three scenarios are handled automatically:
+
+        1. **All penalties zero** ΓÇö only the fitness curve is drawn.
+        2. **No feasible solution found** ΓÇö only the penalty curve is drawn.
+        3. **Mixed** ΓÇö penalty is drawn (shifted) up to the feasibility
+           transition point, then fitness takes over.
+
+        Parameters
+        ----------
+        save_path : str | Path | None
+            If given, the figure is saved at 300 DPI.  Otherwise ``plt.show()``
+            is called.
+        """
+        _apply_academic_style()
+        if not self.iterations_data:
+            print("No data to plot!")
+            return
+
+        first_zero_idx, first_feasible_fitness, all_zero = self._phase_detection()
+        (
+            iterations,
+            fitness_values,
+            penalty_values,
+            min_fitness,
+            last_fitness,
+            max_penalty,
+            first_fitness,
+        ) = self._extract_series()
+
+        (
+            penalties_draw,
+            fitness_draw,
+            min_fitness,
+            last_fitness,
+            y_pen_zero_phys,
+            y_pen_max_phys,
+        ) = self._build_draw_series(
+            all_zero=all_zero,
+            first_zero_idx=first_zero_idx,
+            first_feasible_fitness=first_feasible_fitness,
+            fitness_values=fitness_values,
+            penalty_values=penalty_values,
+            min_fitness=min_fitness,
+            last_fitness=last_fitness,
+            max_penalty=max_penalty,
+        )
+
+        fig, ax1 = plt.subplots(figsize=(12, 6))
+        ax1.set_xlabel(self.x_label)
+        ax2 = ax1.twinx()
 
+        valid_fit_vals = self._plot_fitness_curve(
+            ax1=ax1,
+            iterations=iterations,
+            fitness_draw=fitness_draw,
+        )
+        valid_pen_vals = self._plot_penalty_curve(
+            ax2=ax2,
+            iterations=iterations,
+            penalties_draw=penalties_draw,
+            valid_fit_vals=valid_fit_vals,
+        )
+
+        self._apply_axis_scaling(
+            ax1=ax1,
+            ax2=ax2,
+            valid_fit_vals=valid_fit_vals,
+            valid_pen_vals=valid_pen_vals,
+            min_fitness=min_fitness,
+            first_fitness=first_fitness,
+            all_zero=all_zero,
+            y_pen_zero_phys=y_pen_zero_phys,
+            y_pen_max_phys=y_pen_max_phys,
+        )
+
+        if self.x_lim is not None:
+            ax1.set_xlim(*self.x_lim)
+            ax2.set_xlim(*self.x_lim)
+
+        self._configure_tick_formatters(
+            ax1=ax1,
+            ax2=ax2,
+            first_zero_idx=first_zero_idx,
+            first_feasible_fitness=first_feasible_fitness,
+            all_zero=all_zero,
+            valid_fit_vals=valid_fit_vals,
+        )
+        self._draw_reference_lines(
+            ax1=ax1,
+            ax2=ax2,
+            first_zero_idx=first_zero_idx,
+            all_zero=all_zero,
+            first_feasible_fitness=first_feasible_fitness,
+            min_fitness=min_fitness,
+            last_fitness=last_fitness,
+            y_pen_zero_phys=y_pen_zero_phys,
+        )
+
+        ax1.set_title(self.title, pad=15)
+        self._draw_legend(ax1=ax1, ax2=ax2)
+        self._draw_info_boxes(fig=fig, ax2=ax2)
+
+        fig.tight_layout()
         if save_path:
             plt.savefig(str(save_path), dpi=300, bbox_inches="tight")
             plt.close(fig)
@@ -540,7 +700,7 @@ def plot_multi_run_convergence(
     std = fitness_matrix.std(axis=0)
 
     if ax is None:
-        fig, ax = plt.subplots(figsize=(12, 6))
+        _, ax = plt.subplots(figsize=(12, 6))
         ax.set_title(title, pad=15)
         owns_figure = True
     else:
diff --git a/benchmarks/welded_beam/dependent_space/run.py b/benchmarks/welded_beam/dependent_space/run.py
index 15af873..5d44b90 100644
--- a/benchmarks/welded_beam/dependent_space/run.py
+++ b/benchmarks/welded_beam/dependent_space/run.py
@@ -128,11 +128,11 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
     tau_dash = P / (math.sqrt(2) * x1 * x2)
     tau_dash_dash = M * R / J
     tau = math.sqrt(tau_dash**2 + 2 * tau_dash * tau_dash_dash * (x2 / (2 * R)) + tau_dash_dash**2)
-    Pc = ((4.013 * E * math.sqrt((x3**2 * x4**6) / 36.0)) / L**2) * (1.0 - (x3 / (2 * L)) * math.sqrt(E / (4 * G)))
+    pc = ((4.013 * E * math.sqrt((x3**2 * x4**6) / 36.0)) / L**2) * (1.0 - (x3 / (2 * L)) * math.sqrt(E / (4 * G)))
 
     g1 = tau - TAU_MAX
     g4 = 0.10471 * x1**2 + 0.04811 * x3 * x4 * (14.0 + x2) - 5.0
-    g7 = P - Pc
+    g7 = P - pc
 
     penalty: float = max(0.0, g1) + max(0.0, g4) + max(0.0, g7)
     cost: float = 1.10471 * x1**2 * x2 + 0.04811 * x3 * x4 * (14.0 + x2)
@@ -180,7 +180,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Welded Beam ΓÇö Dependent Space")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.4f}\n" f"Penalty: {result.best_penalty:.4f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.4f}\nPenalty: {result.best_penalty:.4f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/benchmarks/welded_beam/static_penalty/run.py b/benchmarks/welded_beam/static_penalty/run.py
index c341d06..b235d50 100644
--- a/benchmarks/welded_beam/static_penalty/run.py
+++ b/benchmarks/welded_beam/static_penalty/run.py
@@ -101,7 +101,7 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
 
     sigma = 6 * P * L / (x4 * x3**2)
     delta = 4 * P * L**3 / (E * x3**3 * x4)
-    Pc = ((4.013 * E * math.sqrt((x3**2 * x4**6) / 36.0)) / L**2) * (1.0 - (x3 / (2 * L)) * math.sqrt(E / (4 * G)))
+    pc = ((4.013 * E * math.sqrt((x3**2 * x4**6) / 36.0)) / L**2) * (1.0 - (x3 / (2 * L)) * math.sqrt(E / (4 * G)))
 
     # --- All 7 constraints ---
     g1 = tau - TAU_MAX
@@ -110,7 +110,7 @@ def objective(harmony: Dict[str, Any]) -> Tuple[float, float]:
     g4 = 0.10471 * x1**2 + 0.04811 * x3 * x4 * (14.0 + x2) - 5.0
     g5 = 0.125 - x1
     g6 = delta - DELTA_MAX
-    g7 = P - Pc
+    g7 = P - pc
 
     penalty: float = sum(max(0.0, g) for g in [g1, g2, g3, g4, g5, g6, g7])
     cost: float = 1.10471 * x1**2 * x2 + 0.04811 * x3 * x4 * (14.0 + x2)
@@ -158,7 +158,7 @@ def main() -> None:
     plotter = ConvergencePlotter(OUTPUT_DIR / "history_data.csv")
     plotter.set_labels(title="Welded Beam ΓÇö Static Penalty")
     plotter.add_info_box(
-        f"Cost: {result.best_fitness:.4f}\n" f"Penalty: {result.best_penalty:.4f}\n" f"Time: {t_elapsed:.2f}s"
+        f"Cost: {result.best_fitness:.4f}\nPenalty: {result.best_penalty:.4f}\nTime: {t_elapsed:.2f}s"
     )
     plotter.plot(save_path=OUTPUT_DIR / "convergence.png")
 
diff --git a/examples/06_steel_beam_design.py b/examples/06_steel_beam_design.py
index 0d3fca0..bb80da6 100644
--- a/examples/06_steel_beam_design.py
+++ b/examples/06_steel_beam_design.py
@@ -80,8 +80,8 @@ def objective(harmony):
     sec = section_var.decode(idx)
 
     # Convert catalogue values to SI base units
-    wy_m3 = sec.Wy_cm3 * 1e-6  # cm³ → m³
-    iz_m4 = sec.Iy_cm4 * 1e-8  # cm⁴ → m⁴  (strong axis = Iy in catalogue)
+    wy_m3 = sec.wy_cm3 * 1e-6  # cm³ → m³
+    iz_m4 = sec.iy_cm4 * 1e-8  # cm⁴ → m⁴  (strong axis = Iy in catalogue)
     tw_m = sec.tw_mm * 1e-3  # mm  → m
     h_m = sec.h_mm * 1e-3  # mm  → m
 
@@ -135,8 +135,8 @@ if __name__ == "__main__":
     idx = result.best_harmony["section"]
     sec = section_var.decode(idx)
 
-    wy_m3 = sec.Wy_cm3 * 1e-6
-    iz_m4 = sec.Iy_cm4 * 1e-8
+    wy_m3 = sec.wy_cm3 * 1e-6
+    iz_m4 = sec.iy_cm4 * 1e-8
     tw_m = sec.tw_mm * 1e-3
     h_m = sec.h_mm * 1e-3
     fy_kpa = fy * 1e3
@@ -156,4 +156,4 @@ if __name__ == "__main__":
     print(f"  Deflect : ╬┤ = {delta*1000:.1f} mm  Γëñ  L/300 = {L/300*1000:.1f} mm  {'Γ£ô' if delta <= L/300 else 'Γ£ù'}")
     print("\nSection properties:")
     print(f"  h = {sec.h_mm:.0f} mm,  b = {sec.b_mm:.0f} mm")
-    print(f"  A = {sec.A_cm2:.1f} cm²,  Iy = {sec.Iy_cm4:.0f} cm⁴,  Wy = {sec.Wy_cm3:.0f} cm³")
+    print(f"  A = {sec.a_cm2:.1f} cm²,  Iy = {sec.iy_cm4:.0f} cm⁴,  Wy = {sec.wy_cm3:.0f} cm³")
diff --git a/examples/07_rc_section_full.py b/examples/07_rc_section_full.py
index aa59c82..1de4d26 100644
--- a/examples/07_rc_section_full.py
+++ b/examples/07_rc_section_full.py
@@ -88,7 +88,7 @@ space.add("d", Discrete(300.0, 50.0, 900.0))
 rebar_var = ACIRebar(
     d_expr=lambda ctx: ctx["d"] / 1000.0,  # mm → m
     cc_expr=cover,
-    fc=lambda ctx: concrete_var.decode(ctx["grade_idx"]).fck_MPa,
+    fc=lambda ctx: concrete_var.decode(ctx["grade_idx"]).fck_mpa,
     fy=fy,
 )
 space.add("rebar", rebar_var)
@@ -127,7 +127,7 @@ def objective(h):
         return 1e12, 1e6
 
     grade = concrete_var.decode(grade_idx)
-    fc = grade.fck_MPa
+    fc = grade.fck_mpa
     steel_area_mm2 = steel_area(code)
 
     # Moment capacity
@@ -171,12 +171,12 @@ if __name__ == "__main__":
     code = result.best_harmony["rebar"]
 
     grade = concrete_var.decode(grade_idx)
-    fc = grade.fck_MPa
+    fc = grade.fck_mpa
     steel_area_mm2 = steel_area(code)
     phi_mn = moment_capacity(b, d, steel_area_mm2, fc, fy)
 
     print("\nOptimal design:")
-    print(f"  Concrete grade : {grade.name}  (fck = {fc} MPa, Ecm = {grade.Ecm_GPa} GPa)")
+    print(f"  Concrete grade : {grade.name}  (fck = {fc} MPa, Ecm = {grade.ecm_gpa} GPa)")
     print(f"  Section        : b = {b:.0f} mm × d = {d:.0f} mm")
     print(f"  {rebar_var.describe(code)}  (A_s = {steel_area_mm2:.0f} mm┬▓)")
     print("\nStrength check:")
diff --git a/harmonix/logging.py b/harmonix/logging.py
index 50488f4..7078b9b 100644
--- a/harmonix/logging.py
+++ b/harmonix/logging.py
@@ -32,6 +32,7 @@ create them manually.
 from __future__ import annotations
 
 import csv
+import re
 import time
 from collections import OrderedDict
 from pathlib import Path
@@ -61,10 +62,15 @@ def _resolve_path(
     3. ``None`` ΓÇö logging disabled.
     """
     if explicit is not None:
-        return Path(explicit)
+        # Normalize the user-provided path so that subsequent filesystem
+        # operations use a canonical absolute location.
+        return Path(explicit).expanduser().resolve()
     if checkpoint_path is not None:
-        base = Path(checkpoint_path)
-        return base.with_name(base.stem + suffix + ".csv")
+        base = Path(checkpoint_path).expanduser().resolve()
+        # Sanitize the derived stem to avoid accidentally embedding unusual
+        # characters from user-controlled filenames into output paths.
+        safe_stem = re.sub(r"[^a-zA-Z0-9._-]+", "_", base.stem)
+        return base.with_name(f"{safe_stem}{suffix}.csv")
     return None
 
 
diff --git a/harmonix/optimizer.py b/harmonix/optimizer.py
index 3e52b29..a734c21 100644
--- a/harmonix/optimizer.py
+++ b/harmonix/optimizer.py
@@ -41,6 +41,9 @@ Harmony = Dict[str, Any]
 Fitness = float
 Penalty = float
 
+# Hard cap to avoid accidental (or malicious) unbounded iteration counts.
+MAX_ITER_CAP: int = 200_000
+
 
 # ---------------------------------------------------------------------------
 # Result container
@@ -301,6 +304,61 @@ class HarmonySearchOptimizer(ABC):
 
     # --- run setup ---------------------------------------------------------
 
+    def _decide_should_resume(self, *, ckpt: Optional[Path], resume: str) -> bool:
+        if resume == "resume":
+            if ckpt is None or not ckpt.exists():
+                raise FileNotFoundError(f"resume='resume' but checkpoint not found: {ckpt}")
+            return True
+        if resume == "auto":
+            return ckpt is not None and ckpt.exists() and ckpt.stat().st_size > 0
+        if resume == "new":
+            return False
+        raise ValueError(f"resume must be 'auto', 'new', or 'resume'; got {resume!r}.")
+
+    def _wrap_objective(self, *, use_cache: bool, cache_maxsize: int) -> Callable:
+        effective_obj: Callable = self.objective
+        if use_cache:
+            self._cache = EvaluationCache(self.objective, maxsize=cache_maxsize)
+            effective_obj = self._cache
+        else:
+            self._cache = None
+        return effective_obj
+
+    def _resolve_logger_paths(
+        self,
+        *,
+        ckpt: Optional[Path],
+        log_init: bool,
+        init_log_path: Optional[Path],
+        log_evaluations: bool,
+        eval_log_path: Optional[Path],
+        log_history: bool,
+        history_log_path: Optional[Path],
+    ) -> Tuple[Optional[Path], Optional[Path], Optional[Path]]:
+        init_path = _resolve_path(init_log_path, ckpt, "_init") if log_init else None
+        eval_path = _resolve_path(eval_log_path, ckpt, "_evals") if log_evaluations else None
+        hist_path = _resolve_path(history_log_path, ckpt, "_history") if log_history else None
+        return init_path, eval_path, hist_path
+
+    def _init_or_resume_memory(
+        self,
+        *,
+        should_resume: bool,
+        memory_size: int,
+        mode: str,
+        ckpt: Optional[Path],
+        effective_obj: Callable,
+    ) -> int:
+        if should_resume:
+            return self.load_checkpoint(ckpt)  # type: ignore[arg-type]
+
+        self._memory = HarmonyMemory(size=memory_size, mode=mode)
+        for _ in range(memory_size):
+            h = self.space.sample_harmony()
+            f, p = effective_obj(h)
+            self._memory.add(h, float(f), float(p))
+        return 0
+
     def _setup_run(
         self,
         *,
@@ -334,51 +392,25 @@ class HarmonySearchOptimizer(ABC):
             file does not exist.
         """
         ckpt = Path(checkpoint_path) if checkpoint_path else None
-
-        # --- decide whether to resume ---
-        should_resume = False
-        if resume == "resume":
-            if ckpt is None or not ckpt.exists():
-                raise FileNotFoundError(f"resume='resume' but checkpoint not found: {ckpt}")
-            should_resume = True
-        elif resume == "auto":
-            should_resume = ckpt is not None and ckpt.exists() and ckpt.stat().st_size > 0
-        elif resume == "new":
-            should_resume = False
-        else:
-            raise ValueError(f"resume must be 'auto', 'new', or 'resume'; got {resume!r}.")
-
-        # --- initialise or restore memory ---
-        if should_resume:
-            start_iter = self.load_checkpoint(ckpt)
-        else:
-            self._memory = HarmonyMemory(size=memory_size, mode=mode)
-            start_iter = 0
-
-        # --- wrap objective with cache ---
-        effective_obj = self.objective
-        if use_cache:
-            self._cache = EvaluationCache(self.objective, maxsize=cache_maxsize)
-            effective_obj = self._cache
-        else:
-            self._cache = None
-
-        # --- populate memory if starting fresh ---
-        if not should_resume:
-            for _ in range(memory_size):
-                h = self.space.sample_harmony()
-                f, p = effective_obj(h)
-                self._memory.add(h, float(f), float(p))
+        should_resume = self._decide_should_resume(ckpt=ckpt, resume=resume)
+        effective_obj = self._wrap_objective(use_cache=use_cache, cache_maxsize=cache_maxsize)
+        start_iter = self._init_or_resume_memory(
+            should_resume=should_resume,
+            memory_size=memory_size,
+            mode=mode,
+            ckpt=ckpt,
+            effective_obj=effective_obj,
+        )
 
         # --- resolve log paths ---
-        init_path = _resolve_path(Path(init_log_path) if init_log_path else None, ckpt, "_init") if log_init else None
-        eval_path = (
-            _resolve_path(Path(eval_log_path) if eval_log_path else None, ckpt, "_evals") if log_evaluations else None
-        )
-        hist_path = (
-            _resolve_path(Path(history_log_path) if history_log_path else None, ckpt, "_history")
-            if log_history
-            else None
+        init_path, eval_path, hist_path = self._resolve_logger_paths(
+            ckpt=ckpt,
+            log_init=log_init,
+            init_log_path=init_log_path,
+            log_evaluations=log_evaluations,
+            eval_log_path=eval_log_path,
+            log_history=log_history,
+            history_log_path=history_log_path,
         )
 
         logger = RunLogger(
@@ -448,30 +480,7 @@ class Minimization(HarmonySearchOptimizer):
     >>> print(result.best_fitness)
     """
 
-    def optimize(  # type: ignore[override]
-        self,
-        *,
-        memory_size: int = 20,
-        hmcr: float = 0.85,
-        par: float = 0.35,
-        max_iter: int = 5000,
-        bw_max: float = 0.05,
-        bw_min: float = 0.001,
-        resume: str = "auto",
-        checkpoint_path: Optional[Path] = None,
-        checkpoint_every: int = 500,
-        use_cache: bool = False,
-        cache_maxsize: int = 4096,
-        log_init: bool = False,
-        init_log_path: Optional[Path] = None,
-        log_evaluations: bool = False,
-        eval_log_path: Optional[Path] = None,
-        log_history: bool = False,
-        history_log_path: Optional[Path] = None,
-        history_every: int = 1,
-        callback: Optional[Callable[[int, OptimizationResult], None]] = None,
-        verbose: bool = False,
-    ) -> OptimizationResult:
+    def optimize(self, **kwargs: Any) -> OptimizationResult:
         """
         Run the minimisation loop.
 
@@ -515,6 +524,40 @@ class Minimization(HarmonySearchOptimizer):
             Raise :exc:`StopIteration` for early exit.
         verbose : bool
         """
+        memory_size: int = int(kwargs.pop("memory_size", 20))
+        hmcr: float = float(kwargs.pop("hmcr", 0.85))
+        par: float = float(kwargs.pop("par", 0.35))
+        max_iter: int = int(kwargs.pop("max_iter", 5000))
+        bw_max: float = float(kwargs.pop("bw_max", 0.05))
+        bw_min: float = float(kwargs.pop("bw_min", 0.001))
+        resume: str = str(kwargs.pop("resume", "auto"))
+        checkpoint_path_in = kwargs.pop("checkpoint_path", None)
+        checkpoint_path: Optional[Path] = Path(checkpoint_path_in) if checkpoint_path_in is not None else None
+        checkpoint_every: int = int(kwargs.pop("checkpoint_every", 500))
+        use_cache: bool = bool(kwargs.pop("use_cache", False))
+        cache_maxsize: int = int(kwargs.pop("cache_maxsize", 4096))
+        log_init: bool = bool(kwargs.pop("log_init", False))
+        init_log_path_in = kwargs.pop("init_log_path", None)
+        init_log_path: Optional[Path] = Path(init_log_path_in) if init_log_path_in is not None else None
+        log_evaluations: bool = bool(kwargs.pop("log_evaluations", False))
+        eval_log_path_in = kwargs.pop("eval_log_path", None)
+        eval_log_path: Optional[Path] = Path(eval_log_path_in) if eval_log_path_in is not None else None
+        log_history: bool = bool(kwargs.pop("log_history", False))
+        history_log_path_in = kwargs.pop("history_log_path", None)
+        history_log_path: Optional[Path] = Path(history_log_path_in) if history_log_path_in is not None else None
+        history_every: int = int(kwargs.pop("history_every", 1))
+        callback = kwargs.pop("callback", None)
+        verbose: bool = bool(kwargs.pop("verbose", False))
+
+        if kwargs:
+            unexpected = ", ".join(sorted(kwargs))
+            raise TypeError(f"Unexpected keyword argument(s): {unexpected}")
+
+        if max_iter < 0:
+            raise ValueError(f"max_iter must be non-negative; got {max_iter}.")
+        max_iter = min(max_iter, MAX_ITER_CAP)
+        checkpoint_every = max(1, checkpoint_every)
+
         start_iter, logger, effective_obj = self._setup_run(
             memory_size=memory_size,
             mode="min",
@@ -530,6 +573,9 @@ class Minimization(HarmonySearchOptimizer):
             history_log_path=history_log_path,
             history_every=history_every,
         )
+        start_iter = max(0, int(start_iter))
+        if start_iter > max_iter:
+            start_iter = max_iter
         if verbose and start_iter > 0:
             print(f"[HS] Resumed from checkpoint at iteration {start_iter}.")
 
@@ -596,32 +642,38 @@ class Maximization(HarmonySearchOptimizer):
         ``objective(harmony) -> (fitness: float, penalty: float)``
     """
 
-    def optimize(  # type: ignore[override]
-        self,
-        *,
-        memory_size: int = 20,
-        hmcr: float = 0.85,
-        par: float = 0.35,
-        max_iter: int = 5000,
-        bw_max: float = 0.05,
-        bw_min: float = 0.001,
-        resume: str = "auto",
-        checkpoint_path: Optional[Path] = None,
-        checkpoint_every: int = 500,
-        use_cache: bool = False,
-        cache_maxsize: int = 4096,
-        log_init: bool = False,
-        init_log_path: Optional[Path] = None,
-        log_evaluations: bool = False,
-        eval_log_path: Optional[Path] = None,
-        log_history: bool = False,
-        history_log_path: Optional[Path] = None,
-        history_every: int = 1,
-        callback: Optional[Callable[[int, OptimizationResult], None]] = None,
-        verbose: bool = False,
-    ) -> OptimizationResult:
+    def optimize(self, **kwargs: Any) -> OptimizationResult:
         """Run maximisation (see :meth:`Minimization.optimize` for parameter docs)."""
 
+        memory_size: int = int(kwargs.pop("memory_size", 20))
+        hmcr: float = float(kwargs.pop("hmcr", 0.85))
+        par: float = float(kwargs.pop("par", 0.35))
+        max_iter: int = int(kwargs.pop("max_iter", 5000))
+        bw_max: float = float(kwargs.pop("bw_max", 0.05))
+        bw_min: float = float(kwargs.pop("bw_min", 0.001))
+        resume: str = str(kwargs.pop("resume", "auto"))
+        checkpoint_path_in = kwargs.pop("checkpoint_path", None)
+        checkpoint_path: Optional[Path] = Path(checkpoint_path_in) if checkpoint_path_in is not None else None
+        checkpoint_every: int = int(kwargs.pop("checkpoint_every", 500))
+        use_cache: bool = bool(kwargs.pop("use_cache", False))
+        cache_maxsize: int = int(kwargs.pop("cache_maxsize", 4096))
+        log_init: bool = bool(kwargs.pop("log_init", False))
+        init_log_path_in = kwargs.pop("init_log_path", None)
+        init_log_path: Optional[Path] = Path(init_log_path_in) if init_log_path_in is not None else None
+        log_evaluations: bool = bool(kwargs.pop("log_evaluations", False))
+        eval_log_path_in = kwargs.pop("eval_log_path", None)
+        eval_log_path: Optional[Path] = Path(eval_log_path_in) if eval_log_path_in is not None else None
+        log_history: bool = bool(kwargs.pop("log_history", False))
+        history_log_path_in = kwargs.pop("history_log_path", None)
+        history_log_path: Optional[Path] = Path(history_log_path_in) if history_log_path_in is not None else None
+        history_every: int = int(kwargs.pop("history_every", 1))
+        callback = kwargs.pop("callback", None)
+        verbose: bool = bool(kwargs.pop("verbose", False))
+
+        if kwargs:
+            unexpected = ", ".join(sorted(kwargs))
+            raise TypeError(f"Unexpected keyword argument(s): {unexpected}")
+
         wrapped_callback = None
         if callback is not None:
 
@@ -715,31 +767,218 @@ class MultiObjective(HarmonySearchOptimizer):
     ...     print(e.objectives)
     """
 
-    def optimize(  # type: ignore[override]
+    def _mo_init_from_checkpoint(
+        self,
+        *,
+        ckpt: Path,
+        verbose: bool,
+    ) -> Tuple[int, HarmonyMemory, ParetoArchive]:
+        payload = json.loads(ckpt.read_text())
+        start_iter = int(payload["iteration"])
+        memory = HarmonyMemory.from_dict(payload["memory"])
+        archive = ParetoArchive.from_dict(payload["archive"])
+        if verbose:
+            print(f"[MO-HS] Resumed from checkpoint at iteration {start_iter}.")
+        return start_iter, memory, archive
+
+    def _mo_init_fresh(
+        self,
+        *,
+        memory_size: int,
+        archive_size: int,
+        effective_obj: Callable,
+        log_init: bool,
+        logger: RunLogger,
+        ckpt: Optional[Path],
+        verbose: bool,
+    ) -> Tuple[int, HarmonyMemory, ParetoArchive]:
+        memory = HarmonyMemory(size=memory_size, mode="min")
+        archive = ParetoArchive(max_size=archive_size)
+
+        for _ in range(memory_size):
+            h = self.space.sample_harmony()
+            objs, p = effective_obj(h)
+            objs = tuple(float(v) for v in objs)
+            p = float(p)
+            memory.add(h, objs[0], p)
+            if p <= 0:
+                archive.add(h, objs)
+
+        # Log init memory to CSV.
+        if log_init:
+            logger.log_init(
+                harmonies=memory.harmonies,
+                fitnesses=memory._fitness,
+                penalties=memory._penalty,
+            )
+
+        # Write checkpoint after initialisation.
+        if ckpt:
+            payload0 = {
+                "iteration": 0,
+                "memory": memory.to_dict(),
+                "archive": archive.to_dict(),
+            }
+            ckpt.write_text(json.dumps(payload0, indent=2))
+
+        _ = verbose  # reserved for future init logging
+        return 0, memory, archive
+
+    def _mo_init_state(
+        self,
+        *,
+        ckpt: Optional[Path],
+        should_resume: bool,
+        memory_size: int,
+        archive_size: int,
+        effective_obj: Callable,
+        log_init: bool,
+        logger: RunLogger,
+        verbose: bool,
+    ) -> Tuple[int, HarmonyMemory, ParetoArchive]:
+        if should_resume:
+            assert ckpt is not None
+            start_iter, memory, archive = self._mo_init_from_checkpoint(ckpt=ckpt, verbose=verbose)
+        else:
+            start_iter, memory, archive = self._mo_init_fresh(
+                memory_size=memory_size,
+                archive_size=archive_size,
+                effective_obj=effective_obj,
+                log_init=log_init,
+                logger=logger,
+                ckpt=ckpt,
+                verbose=verbose,
+            )
+        self._memory = memory
+        return start_iter, memory, archive
+
+    def _mo_improvise_new_h(
+        self,
+        *,
+        archive: ParetoArchive,
+        hmcr: float,
+        par: float,
+        bw: float,
+    ) -> Harmony:
+        if archive.entries:
+            return self._improvise_from_archive(hmcr, par, archive, bw)
+        return self._improvise(hmcr, par, bw)
+
+    def _mo_evaluate_new_h(
+        self,
+        *,
+        effective_obj: Callable,
+        new_h: Harmony,
+    ) -> Tuple[Tuple[float, ...], float]:
+        objs, p = effective_obj(new_h)
+        objs = tuple(float(v) for v in objs)
+        p = float(p)
+        return objs, p
+
+    def _mo_update_memory_and_archive(
+        self,
+        *,
+        new_h: Harmony,
+        objs: Tuple[float, ...],
+        p: float,
+        archive: ParetoArchive,
+    ) -> None:
+        self._memory.try_replace_worst(new_h, objs[0], p)
+        if p <= 0:
+            archive.add(new_h, objs)
+
+    def _mo_maybe_log_and_callback(
+        self,
+        *,
+        logger: RunLogger,
+        it_plus_1: int,
+        new_h: Harmony,
+        obj0: float,
+        p: float,
+        archive: ParetoArchive,
+        archive_history: List[int],
+        callback: Optional[Callable[[int, ParetoResult], None]],
+        verbose: bool,
+        t0: float,
+    ) -> None:
+        logger.log_evaluation(it_plus_1, new_h, obj0, p)
+
+        if verbose:
+            print(f"[MO-HS] iter {it_plus_1:>6d} | archive = {len(archive):>4d} solutions")
+
+        if callback is not None:
+            partial = ParetoResult(
+                front=archive.front(),
+                archive_history=archive_history,
+                iterations=it_plus_1,
+                elapsed_seconds=time.perf_counter() - t0,
+            )
+            callback(it_plus_1, partial)
+
+    def _mo_maybe_checkpoint(
         self,
         *,
-        memory_size: int = 30,
-        hmcr: float = 0.85,
-        par: float = 0.35,
-        max_iter: int = 5000,
-        bw_max: float = 0.05,
-        bw_min: float = 0.001,
-        archive_size: int = 100,
-        resume: str = "auto",
-        checkpoint_path: Optional[Path] = None,
-        checkpoint_every: int = 500,
-        use_cache: bool = False,
-        cache_maxsize: int = 4096,
-        log_init: bool = False,
-        init_log_path: Optional[Path] = None,
-        log_evaluations: bool = False,
-        eval_log_path: Optional[Path] = None,
-        log_history: bool = False,
-        history_log_path: Optional[Path] = None,
-        history_every: int = 1,
-        callback: Optional[Callable[[int, ParetoResult], None]] = None,
-        verbose: bool = False,
-    ) -> ParetoResult:
+        ckpt: Optional[Path],
+        checkpoint_every: int,
+        it_plus_1: int,
+        archive: ParetoArchive,
+    ) -> None:
+        if ckpt and (it_plus_1 % checkpoint_every == 0):
+            payload_it = {
+                "iteration": it_plus_1,
+                "memory": self._memory.to_dict(),
+                "archive": archive.to_dict(),
+            }
+            ckpt.write_text(json.dumps(payload_it, indent=2))
+
+    def _mo_iteration_update(
+        self,
+        *,
+        it: int,
+        start_iter: int,
+        max_iter: int,
+        hmcr: float,
+        par: float,
+        bw_max: float,
+        bw_min: float,
+        archive: ParetoArchive,
+        effective_obj: Callable,
+        callback: Optional[Callable[[int, ParetoResult], None]],
+        verbose: bool,
+        ckpt: Optional[Path],
+        checkpoint_every: int,
+        logger: RunLogger,
+        archive_history: List[int],
+        t0: float,
+    ) -> None:
+        bw = self._compute_bw(it - start_iter, max_iter - start_iter, bw_max, bw_min)
+        new_h = self._mo_improvise_new_h(archive=archive, hmcr=hmcr, par=par, bw=bw)
+        objs, p = self._mo_evaluate_new_h(effective_obj=effective_obj, new_h=new_h)
+
+        self._mo_update_memory_and_archive(new_h=new_h, objs=objs, p=p, archive=archive)
+        archive_history.append(len(archive))
+
+        it_plus_1 = it + 1
+        self._mo_maybe_log_and_callback(
+            logger=logger,
+            it_plus_1=it_plus_1,
+            new_h=new_h,
+            obj0=objs[0],
+            p=p,
+            archive=archive,
+            archive_history=archive_history,
+            callback=callback,
+            verbose=verbose,
+            t0=t0,
+        )
+        self._mo_maybe_checkpoint(
+            ckpt=ckpt,
+            checkpoint_every=checkpoint_every,
+            it_plus_1=it_plus_1,
+            archive=archive,
+        )
+
+    def optimize(self, **kwargs: Any) -> ParetoResult:
         """
         Run the multi-objective loop.
 
@@ -767,27 +1006,53 @@ class MultiObjective(HarmonySearchOptimizer):
         callback : callable, optional
         verbose : bool
         """
-        # --- setup cache and logger ---
-        effective_obj = self.objective
-        if use_cache:
-            self._cache = EvaluationCache(self.objective, maxsize=cache_maxsize)
-            effective_obj = self._cache
-        else:
-            self._cache = None
-
+        memory_size: int = int(kwargs.pop("memory_size", 30))
+        hmcr: float = float(kwargs.pop("hmcr", 0.85))
+        par: float = float(kwargs.pop("par", 0.35))
+        max_iter: int = int(kwargs.pop("max_iter", 5000))
+        bw_max: float = float(kwargs.pop("bw_max", 0.05))
+        bw_min: float = float(kwargs.pop("bw_min", 0.001))
+        archive_size: int = int(kwargs.pop("archive_size", 100))
+        resume: str = str(kwargs.pop("resume", "auto"))
+        checkpoint_path_in = kwargs.pop("checkpoint_path", None)
+        checkpoint_path: Optional[Path] = Path(checkpoint_path_in) if checkpoint_path_in is not None else None
+        checkpoint_every: int = int(kwargs.pop("checkpoint_every", 500))
+        use_cache: bool = bool(kwargs.pop("use_cache", False))
+        cache_maxsize: int = int(kwargs.pop("cache_maxsize", 4096))
+        log_init: bool = bool(kwargs.pop("log_init", False))
+        init_log_path_in = kwargs.pop("init_log_path", None)
+        init_log_path: Optional[Path] = Path(init_log_path_in) if init_log_path_in is not None else None
+        log_evaluations: bool = bool(kwargs.pop("log_evaluations", False))
+        eval_log_path_in = kwargs.pop("eval_log_path", None)
+        eval_log_path: Optional[Path] = Path(eval_log_path_in) if eval_log_path_in is not None else None
+        log_history: bool = bool(kwargs.pop("log_history", False))
+        history_log_path_in = kwargs.pop("history_log_path", None)
+        history_log_path: Optional[Path] = Path(history_log_path_in) if history_log_path_in is not None else None
+        history_every: int = int(kwargs.pop("history_every", 1))
+        callback = kwargs.pop("callback", None)
+        verbose: bool = bool(kwargs.pop("verbose", False))
+
+        if kwargs:
+            unexpected = ", ".join(sorted(kwargs))
+            raise TypeError(f"Unexpected keyword argument(s): {unexpected}")
+
+        if max_iter < 0:
+            raise ValueError(f"max_iter must be non-negative; got {max_iter}.")
+        max_iter = min(max_iter, MAX_ITER_CAP)
+        checkpoint_every = max(1, checkpoint_every)
+
+        effective_obj = self._wrap_objective(use_cache=use_cache, cache_maxsize=cache_maxsize)
         ckpt = Path(checkpoint_path) if checkpoint_path else None
 
-        # --- resolve log paths ---
-        init_path = _resolve_path(Path(init_log_path) if init_log_path else None, ckpt, "_init") if log_init else None
-        eval_path = (
-            _resolve_path(Path(eval_log_path) if eval_log_path else None, ckpt, "_evals") if log_evaluations else None
-        )
-        hist_path = (
-            _resolve_path(Path(history_log_path) if history_log_path else None, ckpt, "_history")
-            if log_history
-            else None
+        init_path, eval_path, hist_path = self._resolve_logger_paths(
+            ckpt=ckpt,
+            log_init=log_init,
+            init_log_path=init_log_path,
+            log_evaluations=log_evaluations,
+            eval_log_path=eval_log_path,
+            log_history=log_history,
+            history_log_path=history_log_path,
         )
-
         logger = RunLogger(
             variable_names=self.space.names(),
             init_log_path=init_path,
@@ -796,96 +1061,46 @@ class MultiObjective(HarmonySearchOptimizer):
             history_every=history_every,
         )
 
-        # --- decide resume ---
-        archive = ParetoArchive(max_size=archive_size)
-        start_iter = 0
+        should_resume = self._decide_should_resume(ckpt=ckpt, resume=resume)
 
-        should_resume = False
-        if resume == "resume":
-            if ckpt is None or not ckpt.exists():
-                raise FileNotFoundError(f"resume='resume' but checkpoint not found: {ckpt}")
-            should_resume = True
-        elif resume == "auto":
-            should_resume = ckpt is not None and ckpt.exists() and ckpt.stat().st_size > 0
-        elif resume == "new":
-            should_resume = False
-        else:
-            raise ValueError(f"resume must be 'auto', 'new', or 'resume'; got {resume!r}.")
+        start_iter, _, archive = self._mo_init_state(
+            ckpt=ckpt,
+            should_resume=should_resume,
+            memory_size=memory_size,
+            archive_size=archive_size,
+            effective_obj=effective_obj,
+            log_init=log_init,
+            logger=logger,
+            verbose=verbose,
+        )
 
-        if should_resume:
-            payload = json.loads(ckpt.read_text())
-            start_iter = int(payload["iteration"])
-            self._memory = HarmonyMemory.from_dict(payload["memory"])
-            archive = ParetoArchive.from_dict(payload["archive"])
-            if verbose:
-                print(f"[MO-HS] Resumed from checkpoint at iteration {start_iter}.")
-        else:
-            self._memory = HarmonyMemory(size=memory_size, mode="min")
-            for _ in range(memory_size):
-                h = self.space.sample_harmony()
-                objs, p = effective_obj(h)
-                objs = tuple(float(v) for v in objs)
-                p = float(p)
-                self._memory.add(h, objs[0], p)
-                if p <= 0:
-                    archive.add(h, objs)
-
-            if log_init:
-                logger.log_init(
-                    harmonies=self._memory.harmonies,
-                    fitnesses=self._memory._fitness,
-                    penalties=self._memory._penalty,
-                )
-            if ckpt:
-                payload0 = {
-                    "iteration": 0,
-                    "memory": self._memory.to_dict(),
-                    "archive": archive.to_dict(),
-                }
-                ckpt.write_text(json.dumps(payload0, indent=2))
+        start_iter = max(0, int(start_iter))
+        if start_iter > max_iter:
+            start_iter = max_iter
 
         archive_history: List[int] = []
         t0 = time.perf_counter()
 
         try:
             for it in range(start_iter, max_iter):
-                bw = self._compute_bw(it - start_iter, max_iter - start_iter, bw_max, bw_min)
-                if archive.entries:
-                    new_h = self._improvise_from_archive(hmcr, par, archive, bw)
-                else:
-                    new_h = self._improvise(hmcr, par, bw)
-
-                objs, p = effective_obj(new_h)
-                objs = tuple(float(v) for v in objs)
-                p = float(p)
-
-                self._memory.try_replace_worst(new_h, objs[0], p)
-                if p <= 0:
-                    archive.add(new_h, objs)
-
-                archive_history.append(len(archive))
-                logger.log_evaluation(it + 1, new_h, objs[0], p)
-
-                if verbose:
-                    print(f"[MO-HS] iter {it + 1:>6d} | archive = {len(archive):>4d} solutions")
-
-                if callback is not None:
-                    partial = ParetoResult(
-                        front=archive.front(),
-                        archive_history=archive_history,
-                        iterations=it + 1,
-                        elapsed_seconds=time.perf_counter() - t0,
-                    )
-                    callback(it + 1, partial)
-
-                if ckpt and (it + 1) % checkpoint_every == 0:
-                    payload_it = {
-                        "iteration": it + 1,
-                        "memory": self._memory.to_dict(),
-                        "archive": archive.to_dict(),
-                    }
-                    ckpt.write_text(json.dumps(payload_it, indent=2))
-
+                self._mo_iteration_update(
+                    it=it,
+                    start_iter=start_iter,
+                    max_iter=max_iter,
+                    hmcr=hmcr,
+                    par=par,
+                    bw_max=bw_max,
+                    bw_min=bw_min,
+                    archive=archive,
+                    effective_obj=effective_obj,
+                    callback=callback,
+                    verbose=verbose,
+                    ckpt=ckpt,
+                    checkpoint_every=checkpoint_every,
+                    logger=logger,
+                    archive_history=archive_history,
+                    t0=t0,
+                )
         except StopIteration:
             pass
 
diff --git a/harmonix/pareto.py b/harmonix/pareto.py
index 79d4754..14387c3 100644
--- a/harmonix/pareto.py
+++ b/harmonix/pareto.py
@@ -103,7 +103,8 @@ def crowding_distances(
     for m in range(n_obj):
         # Sort by objective m
         # Keep the key lambda single-argument so mypy can infer the type of `i`.
-        order = sorted(range(n), key=lambda i: objective_vectors[i][m])
+        # Capture *m* as a default argument to avoid closure surprises.
+        order = sorted(range(n), key=lambda i, m=m: objective_vectors[i][m])
         f_min = objective_vectors[order[0]][m]
         f_max = objective_vectors[order[-1]][m]
         span = f_max - f_min
diff --git a/harmonix/spaces/engineering.py b/harmonix/spaces/engineering.py
index 28c003f..32678e6 100644
--- a/harmonix/spaces/engineering.py
+++ b/harmonix/spaces/engineering.py
@@ -122,15 +122,23 @@ _DIAMETERS: List[float] = [
 ]
 _COUNTS: List[int] = list(range(4, 42))  # 4 … 41 bars
 
+# All 8 Moore-neighbour offsets for the (diameter × count) discrete grid.
+_MOORE_OFFSETS: List[Tuple[int, int]] = [
+    (di, dj)
+    for di in (-1, 0, 1)
+    for dj in (-1, 0, 1)
+    if not (di == 0 and dj == 0)
+]
+
 
 # ---------------------------------------------------------------------------
 # ACI 318 limit computations  (shared between single- and double-row)
 # ---------------------------------------------------------------------------
 
 
-def _aci_limits(fc: float, fy: float) -> Tuple[float, float, float, float, float]:
+def _aci_limits(fc: float, fy: float) -> Tuple[float, float, float, float]:
     """
-    Return (β₁, φ, ε_cu, ρ_min, ρ_max) using MPa units throughout.
+    Return (β₁, ε_cu, ρ_min, ρ_max) using MPa units throughout.
 
     Parameters
     ----------
@@ -152,7 +160,6 @@ def _aci_limits(fc: float, fy: float) -> Tuple[float, float, float, float, float
     else:
         beta1 = 0.65
 
-    phi = 0.85
     eps_c = 0.003  # ACI crushing strain
 
     rho_min = max(
@@ -161,7 +168,7 @@ def _aci_limits(fc: float, fy: float) -> Tuple[float, float, float, float, float
     )
     rho_max = 0.85 * beta1 * (fc / fy) * (eps_c / (eps_c + 0.004))
 
-    return beta1, phi, eps_c, rho_min, rho_max
+    return beta1, eps_c, rho_min, rho_max
 
 
 def _bar_is_valid_single(
@@ -170,7 +177,6 @@ def _bar_is_valid_single(
     d_eff: float,
     cc: float,
     beta1: float,
-    phi: float,
     eps_c: float,
     rho_min: float,
     rho_max: float,
@@ -269,7 +275,7 @@ class ACIRebar(Variable):
 
     def _valid_codes(self, ctx: Context) -> List[int]:
         d, cc, fc, fy = self._resolve(ctx)
-        beta1, phi, eps_c, rho_min, rho_max = _aci_limits(fc, fy)
+        beta1, eps_c, rho_min, rho_max = _aci_limits(fc, fy)
         codes: List[int] = []
         for i, (area50, dia) in enumerate(zip(_AREAS_50, _DIAMETERS)):
             for j, count in enumerate(_COUNTS):
@@ -279,7 +285,6 @@ class ACIRebar(Variable):
                     d,
                     cc,
                     beta1,
-                    phi,
                     eps_c,
                     rho_min,
                     rho_max,
@@ -313,16 +318,16 @@ class ACIRebar(Variable):
 
         i = value // self._n
         j = value % self._n
+
         neighbours: List[int] = []
-        for di in (-1, 0, 1):
-            for dj in (-1, 0, 1):
-                if di == 0 and dj == 0:
-                    continue
-                ni, nj = i + di, j + dj
-                if 0 <= ni < len(_DIAMETERS) and 0 <= nj < len(_COUNTS):
-                    code = ni * self._n + nj
-                    if code in valid:
-                        neighbours.append(code)
+        for di, dj in _MOORE_OFFSETS:
+            ni, nj = i + di, j + dj
+            if not (0 <= ni < len(_DIAMETERS) and 0 <= nj < len(_COUNTS)):
+                continue
+            code = ni * self._n + nj
+            if code in valid:
+                neighbours.append(code)
+
         return random.choice(neighbours) if neighbours else value
 
     # --- decode ------------------------------------------------------------
@@ -406,7 +411,7 @@ class ACIDoubleRebar(Variable):
 
     def _valid_codes(self, ctx: Context) -> List[int]:
         d1, d2, cc = self._resolve(ctx)
-        beta1, phi, eps_c, rho_min, rho_max = _aci_limits(self.fc, self.fy)
+        beta1, eps_c, rho_min, rho_max = _aci_limits(self.fc, self.fy)
         codes: List[int] = []
         for i, (area50, dia) in enumerate(zip(_AREAS_50, _DIAMETERS)):
             for j, count in enumerate(_COUNTS):
@@ -416,7 +421,6 @@ class ACIDoubleRebar(Variable):
                     d1,
                     cc,
                     beta1,
-                    phi,
                     eps_c,
                     rho_min,
                     rho_max,
@@ -430,7 +434,6 @@ class ACIDoubleRebar(Variable):
                     d2,
                     cc,
                     beta1,
-                    phi,
                     eps_c,
                     rho_min,
                     rho_max,
@@ -456,16 +459,16 @@ class ACIDoubleRebar(Variable):
             return value
         i = value // self._n
         j = value % self._n
+
         neighbours: List[int] = []
-        for di in (-1, 0, 1):
-            for dj in (-1, 0, 1):
-                if di == 0 and dj == 0:
-                    continue
-                ni, nj = i + di, j + dj
-                if 0 <= ni < len(_DIAMETERS) and 0 <= nj < len(_COUNTS):
-                    code = ni * self._n + nj
-                    if code in valid:
-                        neighbours.append(code)
+        for di, dj in _MOORE_OFFSETS:
+            ni, nj = i + di, j + dj
+            if not (0 <= ni < len(_DIAMETERS) and 0 <= nj < len(_COUNTS)):
+                continue
+            code = ni * self._n + nj
+            if code in valid:
+                neighbours.append(code)
+
         return random.choice(neighbours) if neighbours else value
 
     def decode(self, code: int) -> Tuple[float, int]:
@@ -530,27 +533,6 @@ class SectionProperties:
     wz_cm3: float
     mass_kg_m: float
 
-    @property
-    def A_cm2(self) -> float:
-        return self.a_cm2
-
-    @property
-    def Iy_cm4(self) -> float:
-        return self.iy_cm4
-
-    @property
-    def Wy_cm3(self) -> float:
-        return self.wy_cm3
-
-    @property
-    def Iz_cm4(self) -> float:
-        return self.iz_cm4
-
-    @property
-    def Wz_cm3(self) -> float:
-        return self.wz_cm3
-
-
 # ---------------------------------------------------------------------------
 # Built-in catalogue  (EN 10365:2017 + AISC 16th Ed. selection)
 # ---------------------------------------------------------------------------
@@ -811,27 +793,6 @@ class ConcreteGradeProperties:
     ecm_gpa: float
     eps_cu: float
 
-    @property
-    def fck_MPa(self) -> float:
-        return self.fck_mpa
-
-    @property
-    def fck_cube_MPa(self) -> float:
-        return self.fck_cube_mpa
-
-    @property
-    def fcm_MPa(self) -> float:
-        return self.fcm_mpa
-
-    @property
-    def fctm_MPa(self) -> float:
-        return self.fctm_mpa
-
-    @property
-    def Ecm_GPa(self) -> float:
-        return self.ecm_gpa
-
-
 _CONCRETE_GRADES: List[ConcreteGradeProperties] = [
     # name        fck   fck_cube  fcm    fctm   Ecm   eps_cu
     ConcreteGradeProperties("C12/15", 12, 15, 20, 1.57, 27.1, 3.5),
@@ -948,23 +909,6 @@ class SoilProfile:
     dr_pct: float
     vs30_mps: float
 
-    @property
-    def N_lo(self) -> int:
-        return self.n_lo
-
-    @property
-    def N_hi(self) -> int:
-        return self.n_hi
-
-    @property
-    def cu_kPa(self) -> float:
-        return self.cu_kpa
-
-    @property
-    def Dr_pct(self) -> float:
-        return self.dr_pct
-
-
 _SOIL_PROFILES: List[SoilProfile] = [
     # Cohesionless (granular) soils
     SoilProfile("Very loose sand", 0, 4, "ZE", 26, 0, 15, 130),
@@ -1091,11 +1035,6 @@ class SeismicZone:
     SD1: float
     PGA: float
 
-    @property
-    def Ss(self) -> float:
-        return self.ss
-
-
 # Representative TBDY 2018 / AFAD TDTH grid (DD-2, 475-yr return period)
 # ss and S1 values span the range of Turkey's seismic map; site amplification
 # factors Fs and F1 applied for each site class per TBDY 2018 Table 2.3-2.4.
diff --git a/harmonix/variables.py b/harmonix/variables.py
index 2292463..a2f8284 100644
--- a/harmonix/variables.py
+++ b/harmonix/variables.py
@@ -112,9 +112,8 @@ class Continuous(Variable):
 
     def __init__(self, lo, hi):
         # Validate static bounds early (callable bounds checked at sample time)
-        if not callable(lo) and not callable(hi):
-            if float(lo) > float(hi):
-                raise ValueError(f"Continuous: lo ({lo}) must be <= hi ({hi}).")
+        if not callable(lo) and not callable(hi) and float(lo) > float(hi):
+            raise ValueError(f"Continuous: lo ({lo}) must be <= hi ({hi}).")
         self._lo = lo
         self._hi = hi
 
diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py
index 0c55858..be5cd50 100644
--- a/tests/test_edge_cases.py
+++ b/tests/test_edge_cases.py
@@ -75,7 +75,7 @@ class TestVariableEdgeCases:
 
     def test_continuous_lo_eq_hi_filter(self):
         v = Continuous(3.0, 3.0)
-        assert v.filter([2.0, 3.0, 4.0], {}) == [3.0]
+        assert v.filter([2.0, 3.0, 4.0], {}) == pytest.approx([3.0])
 
     def test_continuous_lo_eq_hi_neighbor(self):
         v = Continuous(3.0, 3.0)
@@ -604,8 +604,8 @@ class TestEngineeringPhysics:
         for idx in var._indices:
             props = var.decode(idx)
             assert (
-                abs(props.fcm_MPa - (props.fck_MPa + 8)) < 0.1
-            ), f"{props.name}: fcm={props.fcm_MPa}, fck+8={props.fck_MPa + 8}"
+                abs(props.fcm_mpa - (props.fck_mpa + 8)) < 0.1
+            ), f"{props.name}: fcm={props.fcm_mpa}, fck+8={props.fck_mpa + 8}"
 
     def test_concrete_grade_ecm_formula(self):
         """EC2: Ecm = 22 * (fcm/10)^0.3 GPa."""
@@ -614,10 +614,10 @@ class TestEngineeringPhysics:
         var = ConcreteGrade()
         for idx in var._indices:
             props = var.decode(idx)
-            expected_ecm = 22.0 * (props.fcm_MPa / 10.0) ** 0.3
+            expected_ecm = 22.0 * (props.fcm_mpa / 10.0) ** 0.3
             assert (
-                abs(props.Ecm_GPa - expected_ecm) < 0.5
-            ), f"{props.name}: Ecm={props.Ecm_GPa:.2f}, expected={expected_ecm:.2f}"
+                abs(props.ecm_gpa - expected_ecm) < 0.5
+            ), f"{props.name}: Ecm={props.ecm_gpa:.2f}, expected={expected_ecm:.2f}"
 
     def test_aci_rebar_all_valid_codes_satisfy_rho(self):
         """_valid_codes içindeki her kod _bar_is_valid_single'dan geçmeli."""
@@ -635,7 +635,7 @@ class TestEngineeringPhysics:
         codes = var._valid_codes({})
         assert len(codes) > 0, "Geçerli kod döndürülmedi"
 
-        beta1, phi, eps_c, rho_min, rho_max = _aci_limits(fc, fy)
+        beta1, eps_c, rho_min, rho_max = _aci_limits(fc, fy)
         n_counts = len(_COUNTS)
 
         for code in codes:
@@ -650,7 +650,6 @@ class TestEngineeringPhysics:
                 d_eff,
                 cc,
                 beta1,
-                phi,
                 eps_c,
                 rho_min,
                 rho_max,
@@ -677,6 +676,6 @@ class TestEngineeringPhysics:
         for idx in var._indices:
             sec = var.decode(idx)
             if sec.h_mm > 0:
-                wy_approx = sec.Iy_cm4 / (sec.h_mm / 2 / 10)
-                ratio = wy_approx / sec.Wy_cm3 if sec.Wy_cm3 > 0 else 1.0
-                assert 0.8 <= ratio <= 1.25, f"{sec.name}: Wy={sec.Wy_cm3:.1f}, approx={wy_approx:.1f}"
+                wy_approx = sec.iy_cm4 / (sec.h_mm / 2 / 10)
+                ratio = wy_approx / sec.wy_cm3 if sec.wy_cm3 > 0 else 1.0
+                assert 0.8 <= ratio <= 1.25, f"{sec.name}: Wy={sec.wy_cm3:.1f}, approx={wy_approx:.1f}"
diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py
index 0baafd1..5d3b79e 100644
--- a/tests/test_optimizer.py
+++ b/tests/test_optimizer.py
@@ -41,17 +41,17 @@ class TestHarmonyMemory:
     def test_best_min(self):
         mem = self._make_mem("min")
         h, f, p = mem.best()
-        assert f == 0.0
+        assert math.isclose(f, 0.0, rel_tol=0.0, abs_tol=1e-12)
 
     def test_best_max(self):
         mem = self._make_mem("max")
         h, f, p = mem.best()
-        assert f == 4.0
+        assert math.isclose(f, 4.0, rel_tol=0.0, abs_tol=1e-12)
 
     def test_worst_min(self):
         mem = self._make_mem("min")
         idx = mem.worst_index()
-        assert mem._fitness[idx] == 4.0
+        assert math.isclose(mem._fitness[idx], 4.0, rel_tol=0.0, abs_tol=1e-12)
 
     def test_feasible_replaces_infeasible(self):
         mem = HarmonyMemory(size=3, mode="min")
@@ -239,7 +239,7 @@ class TestMinimization:
             return abs(h["n"] - 7.0), 0.0
 
         result = Minimization(space, obj).optimize(memory_size=10, max_iter=200)
-        assert result.best_harmony["n"] == 7.0
+        assert math.isclose(result.best_harmony["n"], 7.0, rel_tol=0.0, abs_tol=1e-12)
 
     def test_verbose_prints_iteration_progress(self, capsys):
         space = self._sphere_space()
diff --git a/tests/test_spaces.py b/tests/test_spaces.py
index 374b4a7..59d9071 100644
--- a/tests/test_spaces.py
+++ b/tests/test_spaces.py
@@ -125,7 +125,7 @@ class TestNegativeReal:
 
     def test_filter_and_neighbor_clamp(self):
         v = NegativeReal(lo=-10.0)
-        assert v.filter([-10.0, -5.0, 0.0, 1.0], {}) == [-10.0, -5.0]
+        assert v.filter([-10.0, -5.0, 0.0, 1.0], {}) == pytest.approx([-10.0, -5.0])
         for _ in range(20):
             nb = v.neighbor(-1e-9, {})
             assert -10.0 <= nb <= -1e-9
@@ -143,7 +143,7 @@ class TestPositiveReal:
 
     def test_filter_and_neighbor_clamp(self):
         v = PositiveReal(hi=10.0)
-        assert v.filter([-1.0, 0.0, 1e-9, 3.0, 11.0], {}) == [1e-9, 3.0]
+        assert v.filter([-1.0, 0.0, 1e-9, 3.0, 11.0], {}) == pytest.approx([1e-9, 3.0])
         for _ in range(20):
             nb = v.neighbor(10.0, {})
             assert 1e-9 <= nb <= 10.0
@@ -387,7 +387,7 @@ class TestSteelSection:
         idx = var.sample({})
         sec = var.decode(idx)
         assert sec.series == "IPE"
-        assert sec.Iy_cm4 > 0
+        assert sec.iy_cm4 > 0
         assert sec.mass_kg_m > 0
 
     def test_filter(self):
@@ -505,7 +505,7 @@ class TestConcreteGrade:
         for _ in range(30):
             idx = var.sample({})
             g = var.decode(idx)
-            assert 25 <= g.fck_MPa <= 50
+            assert 25 <= g.fck_mpa <= 50
 
     def test_filter(self):
         var = ConcreteGrade()
@@ -527,8 +527,8 @@ class TestConcreteGrade:
         var = ConcreteGrade()
         for idx in var._indices:
             g = var.decode(idx)
-            assert g.fcm_MPa == g.fck_MPa + 8
-            assert g.Ecm_GPa > 0
+            assert g.fcm_mpa == g.fck_mpa + 8
+            assert g.ecm_gpa > 0
 
     def test_invalid_range_raises(self):
         with pytest.raises(ValueError):
diff --git a/tests/test_variables.py b/tests/test_variables.py
index df70b09..de825bd 100644
--- a/tests/test_variables.py
+++ b/tests/test_variables.py
@@ -39,20 +39,20 @@ from harmonix.variables import (
 class TestFrange:
     def test_simple_integer_steps(self):
         result = _frange(0.0, 1.0, 3.0)
-        assert result == [0.0, 1.0, 2.0, 3.0]
+        assert result == pytest.approx([0.0, 1.0, 2.0, 3.0])
 
     def test_endpoint_always_included(self):
         """hi must always appear as the last element."""
         result = _frange(0.0, 3.0, 10.0)
-        assert result[-1] == 10.0
+        assert result[-1] == pytest.approx(10.0)
 
     def test_fractional_step(self):
         result = _frange(0.0, 0.5, 2.0)
-        assert result == [0.0, 0.5, 1.0, 1.5, 2.0]
+        assert result == pytest.approx([0.0, 0.5, 1.0, 1.5, 2.0])
 
     def test_single_element_lo_equals_hi(self):
         result = _frange(5.0, 1.0, 5.0)
-        assert result == [5.0]
+        assert result == pytest.approx([5.0])
 
     def test_step_larger_than_range(self):
         """step > (hi - lo): should return [lo, hi]."""
@@ -77,8 +77,8 @@ class TestFrange:
 
     def test_negative_range(self):
         result = _frange(-3.0, 1.0, 0.0)
-        assert result[0] == -3.0
-        assert result[-1] == 0.0
+        assert result[0] == pytest.approx(-3.0)
+        assert result[-1] == pytest.approx(0.0)
 
 
 # ===========================================================================
@@ -188,7 +188,7 @@ class TestContinuousFilter:
     def test_filter_keeps_in_range(self):
         v = Continuous(2.0, 5.0)
         result = v.filter([1.0, 2.0, 3.5, 5.0, 6.0], {})
-        assert result == [2.0, 3.5, 5.0]
+        assert result == pytest.approx([2.0, 3.5, 5.0])
 
     def test_filter_empty_input(self):
         v = Continuous(0.0, 1.0)
@@ -206,7 +206,7 @@ class TestContinuousFilter:
     def test_filter_with_dependent_bounds(self):
         v = Continuous(lo=lambda ctx: ctx["lo"], hi=lambda ctx: ctx["hi"])
         result = v.filter([0.0, 1.5, 3.0, 5.0], {"lo": 1.0, "hi": 4.0})
-        assert result == [1.5, 3.0]
+        assert result == pytest.approx([1.5, 3.0])
 
     def test_filter_boundary_values_included(self):
         v = Continuous(1.0, 4.0)
@@ -635,7 +635,7 @@ class TestDependentChain:
         b = Continuous(lo=lambda ctx: ctx["a"], hi=10.0)
         ctx = {"a": 5.0}
         result = b.filter([3.0, 5.0, 7.0, 11.0], ctx)
-        assert result == [5.0, 7.0]
+        assert result == pytest.approx([5.0, 7.0])
 
     def test_neighbor_with_dependent_bounds(self):
         """Dependent neighbor [a, 10] aral─▒─ƒ─▒nda kalmal─▒."""
�