"""NSGA-II multi-objective Pareto optimisation of per-element c-factors.
Replaces the exhaustive ND grid scan (``HiFiWorkflow.scan_cfactor_grid``) as
the default multi-element c-optimisation path. The grid scan remains practical
for N_el ≤ 2.
Public API
----------
HiFiCFactorProblem — batch pymoo Problem evaluating scan-axis vectors
optimize_c_factors_nsga2 — run NSGA-II, return result + history + HV trace
fidelity_envelope — maximum fidelity at each subset size
extract_pareto_front — post-process result into DataFrame + auto-pick point
plot_nsga2_results — produce NSGA-II diagnostic plots
"""
from __future__ import annotations
import logging
import warnings
from pathlib import Path
from typing import Dict, List, Optional, Sequence
logger = logging.getLogger(__name__)
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable, viridis as _cm_viridis
from matplotlib.colors import Normalize
import numpy as np
import pandas as pd
from scipy.signal import find_peaks
from scipy.stats import qmc
import contextlib
from joblib import Parallel, delayed
from pymoo.algorithms.moo.nsga2 import NSGA2
from pymoo.core.callback import Callback
from pymoo.core.duplicate import DefaultDuplicateElimination
from pymoo.core.problem import Problem
from pymoo.indicators.hv import HV
from pymoo.optimize import minimize
from .hifi_workflow import (
_count_kind_k_from_fraction,
_fit_cfactor_worker,
_shm_scope,
_strip_lib_for_workers,
)
from .cluster_backends import BACKEND_REGISTRY
# ─────────────────────────────────────────────────────────────────────────────
# pymoo Problem (batch, parallel)
# ─────────────────────────────────────────────────────────────────────────────
[docs]
class HiFiCFactorProblem(Problem):
"""Batch pymoo Problem wrapping _fit_cfactor_worker with joblib parallelism.
Evaluates each generation's points that miss the cross-generation cache
(count-kind backends also skip in-batch duplicates) serially
(``n_jobs == 1``) or in one ``joblib.Parallel`` call; workers get a
symbols-only library when ``n_jobs > 1``. Descriptors and PCA embeddings
come from SharedMemory when ``optimize_c_factors_nsga2`` runs with its own
``n_jobs > 1`` (the pattern of ``HiFiWorkflow._scan_cfactor``).
Objectives (both minimised):
F[0] = -L_min (maximise likeness)
F[1] = R (minimise redundancy)
"""
def __init__(
self,
species_order: List[str],
xl: Sequence[float],
xu: Sequence[float],
evaluator_kwargs: dict,
n_jobs: int = 1,
):
super().__init__(
n_var=len(species_order),
n_obj=2,
n_ieq_constr=0,
xl=np.array(xl, dtype=float),
xu=np.array(xu, dtype=float),
)
self.species_order = species_order
self.evaluator_kwargs = evaluator_kwargs
self.n_jobs = int(n_jobs)
self.cluster_backend = evaluator_kwargs.get('cluster_backend', 'decaf')
# Same registry lookup the worker body uses for f→k conversion; drives
# the cache-key granularity (count-kind keys on realized k, not f).
_caps = (BACKEND_REGISTRY[self.cluster_backend]
if isinstance(self.cluster_backend, str) else self.cluster_backend)
self._granularity_kind = _caps.granularity_kind
self._cache: dict = {}
self.history: list = []
self._shm_meta = None # injected by optimize_c_factors_nsga2 before minimize
lib = evaluator_kwargs['lib']
lib_workers = _strip_lib_for_workers(lib) if self.n_jobs > 1 else lib
# smap_per_el → struct_map_per_el rename for _fit_cfactor_worker
self._worker_kwargs = dict(
species = species_order,
lib = lib_workers,
lib_path = evaluator_kwargs.get('lib_path', ''),
desc_per_el = evaluator_kwargs['desc_per_el'],
struct_map_per_el = evaluator_kwargs['smap_per_el'],
mu_per_el = evaluator_kwargs['mu_per_el'],
counts_L_per_el = evaluator_kwargs['counts_L_per_el'],
k_occ_L_per_el = evaluator_kwargs['k_occ_L_per_el'],
H_L_per_el = evaluator_kwargs['H_L_per_el'],
eps_noise_per_el = evaluator_kwargs['eps_noise_per_el'],
embed_op_per_el = evaluator_kwargs['embed_op_per_el'],
cluster_backend = self.cluster_backend,
min_cluster_size = evaluator_kwargs.get('min_cluster_size', 2),
)
[docs]
def canonical_key(self, x):
"""Canonical cache key for a scan-axis vector x.
Count-kind: the realized k-vector (round(f·k_micro) per element, same
conversion the worker performs internally) — many f-vectors map to one
backend configuration. Distance-kind: round-to-4dp f-vector (legacy).
"""
if self._granularity_kind == 'count':
return tuple(
_count_kind_k_from_fraction(
v, len(self._worker_kwargs['counts_L_per_el'][el]))
for v, el in zip(x, self.species_order)
)
return tuple(round(float(v), 4) for v in x)
def _evaluate(self, X, out, *args, **kwargs):
n = X.shape[0]
penalty = [0.0, 1.0] # worst-case: L_min=0, R=1 (maximum redundancy)
# Partition into cached / uncached. Canonical key, with an in-batch
# `seen` set for count-kind only: many count-kind f-vectors realize to
# the same k-vector within one generation, so skip duplicates before
# worker dispatch. Distance-kind keeps its pre-existing behavior
# exactly (each population member gets its own worker call unless it
# hits the cross-generation self._cache, as before) — an in-batch
# dedup there would reuse one near-duplicate's result for another
# point that was previously evaluated independently (issue #82
# codex-review finding: this is a real behavior change, not covered
# by the "additive for count-kind only" design promise).
dedup_in_batch = self._granularity_kind == 'count'
uncached_idxs: List[int] = []
uncached_c_dicts: List[dict] = []
seen: set = set()
for i, x in enumerate(X):
key = self.canonical_key(x)
if key in self._cache or (dedup_in_batch and key in seen):
continue
if dedup_in_batch:
seen.add(key)
uncached_idxs.append(i)
uncached_c_dicts.append(
{el: float(x[j]) for j, el in enumerate(self.species_order)}
)
# Evaluate uncached points (serial or parallel)
if uncached_c_dicts:
if self._shm_meta is not None:
meta_desc, meta_emb = self._shm_meta
worker_kwargs = {
**self._worker_kwargs,
'desc_per_el': None,
'_shm_meta_desc': meta_desc,
'_shm_meta_emb': meta_emb,
}
else:
worker_kwargs = self._worker_kwargs
if self.n_jobs == 1:
raw_rows = [_fit_cfactor_worker(c, **worker_kwargs)
for c in uncached_c_dicts]
else:
raw_rows = Parallel(n_jobs=self.n_jobs, prefer='processes')(
delayed(_fit_cfactor_worker)(
c, **worker_kwargs, worker_subprocess=True)
for c in uncached_c_dicts
)
# Augment + cache + record history
for i_unc, row in zip(uncached_idxs, raw_rows):
if row is None:
continue
# Scan workers expose fitted per-element cache payloads for the
# line/grid scan. NSGA-II never consumes them; retaining one
# O(N) payload per evaluated vector makes history scale with
# the full corpus instead of scalar objectives.
row.pop('_el_fits', None)
x = X[i_unc]
row['c_vec'] = np.array(x, dtype=float)
row['delta_vec'] = np.array(
[row[f'delta_{el}'] for el in self.species_order],
dtype=float,
)
row['L_min'] = row['L']
row['patient_OOD'] = 0
key = self.canonical_key(x)
row['_cache_key'] = key
self._cache[key] = row
self.history.append(row)
# Assemble out['F'] (pymoo API: objectives go in out['F'])
F_out = np.empty((n, 2))
for i, x in enumerate(X):
key = self.canonical_key(x)
result = self._cache.get(key)
F_out[i] = (
[-result['L_min'], result['R']]
if result is not None and np.isfinite(result.get('R', float('nan')))
else penalty
)
out['F'] = F_out
# ─────────────────────────────────────────────────────────────────────────────
# HV convergence callback
# ─────────────────────────────────────────────────────────────────────────────
class _HVCallback(Callback):
def __init__(self, ref_point: np.ndarray, bar=None):
super().__init__()
self.ref_point = ref_point
self.hv_history: List[float] = []
self.bar = bar # tqdm generation bar (or None)
def notify(self, algorithm):
F = algorithm.opt.get('F')
if F is not None and len(F):
self.hv_history.append(float(HV(ref_point=self.ref_point).do(F)))
if self.bar is not None and self.bar.n < (self.bar.total or float('inf')):
self.bar.update(1)
# ─────────────────────────────────────────────────────────────────────────────
# Optimizer driver
# ─────────────────────────────────────────────────────────────────────────────
[docs]
def optimize_c_factors_nsga2(
problem: HiFiCFactorProblem,
n_gen: int = 50,
pop_size: int = 40,
seed: int = 42,
c_star_1d: Optional[float] = None,
c_vec_2d_best: Optional[np.ndarray] = None,
n_jobs: int = 1,
):
"""Run NSGA-II on *problem* and return results.
Parameters
----------
problem : HiFiCFactorProblem
n_gen, pop_size, seed : NSGA-II budget parameters.
Rule of thumb: pop_size ≥ 10 × N_focus_elements (min 40).
Watch pareto_convergence.png — increase n_gen if HV hasn't plateaued.
c_star_1d : float or None
Common 1-D c* from the 1-D scan; seeds the initial population.
c_vec_2d_best : array-like or None
Best c-vector from a prior 2-D grid scan (if available).
n_jobs : int
When > 1, descriptors and their PCA embeddings are copied to
SharedMemory once for the whole run and the problem's direct
descriptor references are dropped before minimisation. Evaluation
parallelism is ``problem.n_jobs``. On every path the direct
descriptor references are dropped after minimisation completes.
Returns
-------
(res, history, hv_history)
res — pymoo Result object (.X, .F on the Pareto front)
history — the problem's list of result dicts, one per worker call
that returned a row (cross-generation cache hits and count-kind
in-batch duplicates are not re-evaluated; failed workers add none)
hv_history — list of hypervolume values per generation
"""
xl = problem.xl
xu = problem.xu
n_var = problem.n_var
# ── Seeded starting points ────────────────────────────────────────────────
seeded: List[np.ndarray] = []
if c_star_1d is not None:
seeded.append(np.full(n_var, float(c_star_1d)))
if c_vec_2d_best is not None:
v = np.asarray(c_vec_2d_best, dtype=float)
if v.shape == (n_var,):
seeded.append(np.clip(v, xl, xu))
# Corners
seeded.append(xl.copy())
seeded.append(xu.copy())
# Diagonal sweep: 8 evenly-spaced points along xl→xu
for t in np.linspace(0.0, 1.0, 8):
seeded.append(xl + t * (xu - xl))
# Deduplicate seeded points (rounded to 4dp)
seen_keys: set = set()
unique_seeded: List[np.ndarray] = []
for pt in seeded:
k = tuple(round(float(v), 4) for v in pt)
if k not in seen_keys:
seen_keys.add(k)
unique_seeded.append(pt)
# LHS fill to reach pop_size. Degenerate dims (xl[i] == xu[i], e.g. a
# small-k_micro count-kind element pinned to floor) are skipped — qmc.scale
# raises on zero-width bounds — and stay at xl.
n_seeded = len(unique_seeded)
n_lhs = max(0, pop_size - n_seeded)
if n_lhs > 0:
free = xu > xl
X_lhs = np.repeat(xl[None, :], n_lhs, axis=0)
if free.any():
sampler = qmc.LatinHypercube(d=int(free.sum()), seed=seed)
X_lhs[:, free] = qmc.scale(sampler.random(n_lhs), xl[free], xu[free])
X_init = np.vstack([np.array(unique_seeded), X_lhs])
else:
X_init = np.array(unique_seeded[:pop_size])
X_init = np.clip(X_init, xl, xu)
# ── Reference point for HV ────────────────────────────────────────────────
# F[0] = -L_min ∈ [-1, 0]; worst case = 0.0 (L_min = 0).
# F[1] = R, treated as a fraction; failed evaluations get penalty R = 1.0.
# pymoo minimises; a front point adds hypervolume only where it is better
# than this reference point in both objectives. R sits 0.01 above the
# failure penalty R = 1.0 (`penalty` in HiFiCFactorProblem._evaluate), so
# real points with L_min > 0 and R ≤ 1.0 count; change the two together.
ref_point = np.array([0.0, 1.0 + 0.01])
from tqdm.auto import tqdm
from atom_hifi._progress import bar_disabled, bar_file
bar = tqdm(total=n_gen, desc='NSGA-II', file=bar_file(),
disable=bar_disabled())
cb = _HVCallback(ref_point, bar)
logger.info(f'\nNSGA-II: N_el={n_var} pop_size={pop_size} n_gen={n_gen}'
f' → ~{n_gen * pop_size} evaluations (seeded={n_seeded}, LHS={n_lhs})'
f' n_jobs={n_jobs}')
# Open shared memory ONCE before minimize so it persists across all generations.
# Workers in each generation's Parallel call attach to the same shm blocks.
use_shm = n_jobs > 1
cm = (_shm_scope(problem.evaluator_kwargs['desc_per_el'], problem.species_order)
if use_shm else contextlib.nullcontext((None, None)))
try:
with cm as (meta_desc, meta_emb):
if use_shm and meta_desc is not None:
problem._shm_meta = (meta_desc, meta_emb)
# From this point every evaluation attaches through metadata.
# Keeping either problem dictionary populated would duplicate
# the descriptor corpus in the main process for all generations.
problem.evaluator_kwargs['desc_per_el'] = None
problem._worker_kwargs['desc_per_el'] = None
algo = NSGA2(
pop_size=pop_size,
sampling=X_init,
eliminate_duplicates=DefaultDuplicateElimination(epsilon=1e-4),
)
res = minimize(
problem, algo,
('n_gen', n_gen),
seed=seed,
callback=cb,
verbose=logger.isEnabledFor(logging.DEBUG),
)
# pymoo's Result retains Problem through both res.problem and
# res.algorithm.problem. Evaluations are complete, so direct
# descriptors have reached their last consumer even on the serial
# path (parallel cleared them immediately after the SHM copy).
problem.evaluator_kwargs['desc_per_el'] = None
problem._worker_kwargs['desc_per_el'] = None
finally:
problem._shm_meta = None
bar.close()
return res, problem.history, cb.hv_history
# ─────────────────────────────────────────────────────────────────────────────
# Post-processing
# ─────────────────────────────────────────────────────────────────────────────
[docs]
def fidelity_envelope(df_pareto: pd.DataFrame,
x_col: str = 'n_structs_S') -> pd.DataFrame:
"""Upper envelope of fidelity F over x_col (default: ``|S|``).
Groups by x_col, keeps max-F row per group, returns sorted ascending.
Rows with non-finite F are dropped; returns empty DataFrame if none remain.
"""
if 'F' not in df_pareto.columns or x_col not in df_pareto.columns:
return df_pareto.iloc[0:0].copy()
df = df_pareto[np.isfinite(df_pareto['F']) & np.isfinite(df_pareto[x_col])]
if df.empty:
return df
idx = df.groupby(x_col)['F'].idxmax()
return df.loc[idx].sort_values(x_col)
# ─────────────────────────────────────────────────────────────────────────────
# Plots
# ─────────────────────────────────────────────────────────────────────────────
def _plot_fidelity_vs_x(
df_pareto, hdf, fb, sel,
x_col, x_label, x_short,
title, l_tol, species_order, save_path,
):
"""Render F vs log(x_col) operating-point plot and save to save_path.
Shared by the |S| (structure-count) and atoms-in-S variants produced by
plot_nsga2_results. ``hdf``, ``fb``, ``sel`` are pre-computed by the
caller so both plots share the same operating-point extraction.
"""
if x_col not in df_pareto.columns:
return
df_x = df_pareto[
np.isfinite(df_pareto['F'])
& np.isfinite(df_pareto[x_col])
& (df_pareto[x_col] > 0)
]
if df_x.empty:
return
def _valid_operating_point(row):
return (
row is not None
and x_col in row.index
and np.isfinite(row[x_col])
and row[x_col] > 0
)
fb_x = fb if _valid_operating_point(fb) else None
sel_x = sel if _valid_operating_point(sel) else None
fig, ax = plt.subplots(figsize=(8, 5.5))
if not hdf.empty and x_col in hdf.columns:
hdf_x = hdf[
np.isfinite(hdf['F'])
& np.isfinite(hdf[x_col])
& (hdf[x_col] > 0)
]
if not hdf_x.empty:
ax.scatter(hdf_x[x_col], hdf_x['F'], c='#dddddd', s=8, zorder=0,
label=f'all history (n={len(hdf_x)})')
h_env = fidelity_envelope(hdf_x, x_col=x_col)
ax.plot(h_env[x_col], h_env['F'], '-', color='#888888',
lw=1.0, zorder=2, label='envelope: all points')
ax.scatter(df_x[x_col], df_x['F'], c='#999999', s=14,
zorder=2, label='front points')
env = fidelity_envelope(df_x, x_col=x_col)
ax.plot(env[x_col], env['F'], 'o-', color='steelblue', zorder=3,
label='F envelope (front)')
env_s = env.sort_values(x_col).reset_index(drop=True)
F_arr = env_s['F'].values
glob_i = int(np.argmax(F_arr))
if len(F_arr) >= 3:
prom = max(0.02, 0.05 * float(np.ptp(F_arr)))
peak_idx, props = find_peaks(F_arr, prominence=prom)
avoid_X = [float(r[x_col]) for r in (fb_x, sel_x) if r is not None]
c_min = {el: df_x[f'c_{el}'].min() for el in species_order}
c_rng = {el: max(df_x[f'c_{el}'].max() - c_min[el], 1e-12)
for el in species_order}
kept = []
for pi in peak_idx[np.argsort(props['prominences'])[::-1]]:
if pi == glob_i:
continue
sp = float(env_s.iloc[pi][x_col])
if any(0.87 < sp / a < 1.15 for a in avoid_X):
continue
kept.append(int(pi))
if len(kept) >= 4:
break
for rank, pi in enumerate(sorted(kept)):
prow = env_s.iloc[pi]
top_el = max(species_order,
key=lambda el: (prow[f'c_{el}'] - c_min[el]) / c_rng[el])
ax.scatter([prow[x_col]], [prow['F']], marker='o', s=70,
facecolors='none', edgecolors='darkorange',
linewidths=1.8, zorder=4)
above = rank % 2 == 0
ax.annotate(
f'{x_short}={int(prow[x_col])}, F={prow["F"]:.3f}, ↑{top_el}',
(prow[x_col], prow['F']), fontsize=7,
ha='center', va='bottom' if above else 'top',
color='darkorange',
xytext=(0, 9 if above else -9), textcoords='offset points')
if fb_x is not None:
ax.scatter([fb_x[x_col]], [fb_x['F']], marker='*', s=240,
color='gold', edgecolors='black', linewidths=0.5, zorder=6,
label=f'fidelity_best (F={fb_x["F"]:.3f}, {x_short}={int(fb_x[x_col])})')
if sel_x is not None:
_sel_tag = f'L_min≥{l_tol}' if l_tol is not None else 'fidelity_best'
ax.scatter([sel_x[x_col]], [sel_x['F']], marker='v', s=150,
color='crimson', edgecolors='black', linewidths=0.5, zorder=5,
label=f'auto-pick ({_sel_tag}, {x_short}={int(sel_x[x_col])})')
ax.set_xscale('log')
ax.set_xlabel(x_label)
ax.set_ylabel('F = L / R')
ax.set_title(title)
ax.legend(fontsize=8)
fig.tight_layout()
fig.savefig(save_path, dpi=150)
plt.close(fig)
[docs]
def plot_nsga2_results(
df_pareto: pd.DataFrame,
history: list,
hv_history: List[float],
species_order: List[str],
best_c_per_el: dict,
l_tol: float = 0.90,
save_dir=None,
):
"""Produce NSGA-II diagnostic plots and save to save_dir."""
save_dir = Path(save_dir) if save_dir is not None else Path('.')
if df_pareto.empty:
logger.warning('plot_nsga2_results: df_pareto is empty — skipping plots.')
return
# Build history arrays for the gray scatter.
hist_R = [e.get('R', float('nan')) for e in history if e is not None]
hist_L = [e.get('L_min', float('nan')) for e in history if e is not None]
# ── Plot A: Pareto front (L vs R) ─────────────────────────────────────────
fig, ax = plt.subplots(figsize=(7, 5))
# Gray cloud — all evaluated points.
_hr = np.array(hist_R, dtype=float)
_hl = np.array(hist_L, dtype=float)
_ok = np.isfinite(_hr) & np.isfinite(_hl)
if _ok.any():
ax.scatter(_hr[_ok], _hl[_ok], c='#cccccc', s=10, zorder=1,
label='all evaluated')
# Pareto front — df_pareto is non-dominated in (L, R) space because the
# NSGA-II objectives are (-L_min, R). Sorting by R gives a monotone curve
# with no loops; no manual filtering needed.
_fr = df_pareto[np.isfinite(df_pareto['R'])].sort_values('R').reset_index(drop=True)
ax.plot(_fr['R'].values, _fr['L_min'].values, 'o-', color='steelblue',
zorder=3, label='Pareto front')
# Annotate a thinned subset with |S|.
_step = max(1, int(np.ceil(len(_fr) / 8)))
for _i, _arow in _fr.iterrows():
if _i % _step != 0 and _i != len(_fr) - 1:
continue
ax.annotate(f'|S|={int(_arow["n_structs_S"])}',
(_arow['R'], _arow['L_min']),
fontsize=7, ha='left', va='bottom',
xytext=(4, 4), textcoords='offset points')
# Gold star = fidelity_best (global argmax F); red triangle = auto-pick
# (smallest R, L_min≥l_tol) — markers consistent with pareto_fidelity.png.
fb = (df_pareto[df_pareto['fidelity_best']].iloc[0]
if 'fidelity_best' in df_pareto.columns and df_pareto['fidelity_best'].any()
else None)
if fb is not None and np.isfinite(fb['R']):
ax.scatter([fb['R']], [fb['L_min']], marker='*', s=240, color='gold',
edgecolors='black', linewidths=0.5, zorder=6,
label=f'fidelity_best (|S|={int(fb["n_structs_S"])})')
sel = (df_pareto[df_pareto['selected']].iloc[0]
if 'selected' in df_pareto.columns and df_pareto['selected'].any()
else None)
pick = sel
if pick is None:
_m = df_pareto[np.all([np.isclose(df_pareto[f'c_{el}'], best_c_per_el[el], atol=1e-4)
for el in species_order], axis=0)]
pick = _m.iloc[0] if not _m.empty else None
if pick is not None and np.isfinite(pick['R']):
ax.scatter([pick['R']], [pick['L_min']], marker='v', s=150, color='crimson',
edgecolors='black', linewidths=0.5, zorder=5,
label=f'auto-pick (|S|={int(pick["n_structs_S"])})')
if l_tol is not None:
ax.axhline(l_tol, color='red', linestyle='--', linewidth=0.8,
label=f'L = {l_tol}')
ax.set_xlabel('R (redundancy)')
ax.set_ylabel('L (likeness)')
ax.set_title('NSGA-II Pareto front (L vs R)')
ax.legend(fontsize=8)
fig.tight_layout()
fig.savefig(save_dir / 'pareto_front.png', dpi=150)
plt.close(fig)
# ── Plots D/E: Fidelity F vs log|S| and F vs log(atoms in S) ────────────
hdf = pd.DataFrame([e for e in history if e is not None])
if not hdf.empty and 'F' in hdf.columns:
hdf = hdf[np.isfinite(hdf['F'])]
else:
hdf = pd.DataFrame()
_plot_fidelity_vs_x(
df_pareto, hdf, fb, sel,
x_col='n_structs_S', x_label='|S| (structures)', x_short='|S|',
title='Fidelity vs log|S| (Pareto front)',
l_tol=l_tol, species_order=species_order,
save_path=save_dir / 'pareto_fidelity.png',
)
if 'n_atoms_S' in df_pareto.columns:
_plot_fidelity_vs_x(
df_pareto, hdf, fb, sel,
x_col='n_atoms_S', x_label='atoms in S', x_short='atoms',
title='Fidelity vs log(atoms in S) (Pareto front)',
l_tol=l_tol, species_order=species_order,
save_path=save_dir / 'pareto_fidelity_atoms.png',
)
# ── Plot B: Parallel coordinates (N_el ≥ 2) ──────────────────────────────
# Standard single-axes parallel coordinates: one labelled vertical axis per
# variable, each normalised to its own min→max (printed at top/bottom), lines
# coloured by L_min. Much clearer than one cramped subplot per axis-pair.
if len(species_order) >= 2:
cols = [f'c_{el}' for el in species_order] + ['L', 'R', 'F', 'n_structs_S']
K = len(cols)
data = df_pareto[cols].to_numpy(dtype=float)
cmin, cmax = data.min(axis=0), data.max(axis=0)
crange = np.where(cmax > cmin, cmax - cmin, 1.0)
norm_data = (data - cmin) / crange # each column → [0, 1]
fmin_vals = df_pareto['L_min'].values
cmap = _cm_viridis
norm = Normalize(fmin_vals.min(), fmin_vals.max())
fig, ax = plt.subplots(figsize=(1.7 * K + 1.5, 5))
xs = np.arange(K)
for r in range(norm_data.shape[0]):
ax.plot(xs, norm_data[r], color=cmap(norm(fmin_vals[r])),
alpha=0.6, lw=1.2)
for j in range(K):
ax.axvline(j, color='0.85', lw=1, zorder=0)
ax.text(j, 1.03, f'{cmax[j]:.3g}', ha='center', va='bottom', fontsize=8)
ax.text(j, -0.03, f'{cmin[j]:.3g}', ha='center', va='top', fontsize=8)
ax.set_xticks(xs)
ax.set_xticklabels(cols, fontsize=10)
ax.set_yticks([])
ax.set_ylim(-0.10, 1.12)
ax.set_xlim(-0.3, K - 0.7)
ax.set_title('Pareto front — parallel coordinates (each axis min→max)')
fig.colorbar(ScalarMappable(cmap=cmap, norm=norm),
ax=ax, label='L_min', pad=0.01)
fig.tight_layout()
fig.savefig(save_dir / 'pareto_parallel.png', dpi=150)
plt.close(fig)
# ── Plot C: HV convergence ────────────────────────────────────────────────
if hv_history:
fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(range(1, len(hv_history) + 1), hv_history, color='steelblue')
ax.set_xlabel('Generation')
ax.set_ylabel('Hypervolume')
ax.set_title('NSGA-II convergence')
fig.tight_layout()
fig.savefig(save_dir / 'pareto_convergence.png', dpi=150)
plt.close(fig)
# ─────────────────────────────────────────────────────────────────────────────
# Sanity summary printer
# ─────────────────────────────────────────────────────────────────────────────
def _print_sanity_summary(
res,
df_pareto: pd.DataFrame,
hv_history: List[float],
species_order: List[str],
best_c_per_el: dict,
l_tol: float,
):
if df_pareto.empty:
logger.warning('_print_sanity_summary: df_pareto is empty — no front to summarise.')
return
n_evals = (res.algorithm.evaluator.n_eval
if hasattr(res, 'algorithm') and res.algorithm is not None else '?')
n_front = len(df_pareto)
hv_final = hv_history[-1] if hv_history else float('nan')
hv_delta = (hv_history[-1] - hv_history[-11]
if len(hv_history) > 10 else float('nan'))
best_fmin_row = df_pareto.loc[df_pareto['L_min'].idxmax()]
best_s_row = df_pareto.iloc[0]
qualifying = (df_pareto[df_pareto['L_min'] >= l_tol]
if l_tol is not None else df_pareto.iloc[0:0])
f95_row = qualifying.iloc[0] if not qualifying.empty else None
def _fmt_row(row):
c_str = '[' + ', '.join(f'{row[f"c_{el}"]:.4f}' for el in species_order) + ']'
return f'c={c_str} L_min={row["L_min"]:.4f} |S|={int(row["n_structs_S"])}'
logger.info('\nOptimizer summary:')
logger.info(f' Total evaluations : {n_evals}')
logger.info(f' Front size : {n_front}')
logger.info(f' Hypervolume : {hv_final:.6f} '
f'(delta over last 10 gens: {hv_delta:+.6f})')
logger.info(f'\n Best L_min point : {_fmt_row(best_fmin_row)}')
if l_tol is None:
logger.info(' L floor disabled (l_tol: null) — auto-pick = fidelity_best')
elif f95_row is not None:
logger.info(f' F ≥ {l_tol} constrained : {_fmt_row(f95_row)}')
else:
logger.info(f' F ≥ {l_tol} constrained : NONE — see warning above')
logger.info(f' Smallest R on front : {_fmt_row(best_s_row)}')
if 'fidelity_best' in df_pareto.columns and df_pareto['fidelity_best'].any():
fb = df_pareto[df_pareto['fidelity_best']].iloc[0]
logger.info(f' Best fidelity (global): {_fmt_row(fb)} F={fb["F"]:.4f}')
if f95_row is not None:
logger.info(f'\n L={f95_row["L"]:.4f} R={f95_row["R"]:.4f} F={f95_row["F"]:.4f} '
f'at constrained optimum')
if l_tol is not None and best_fmin_row['L_min'] < l_tol:
logger.warning(f'best L_min={best_fmin_row["L_min"]:.4f} < {l_tol}. '
'Check eps_noise calibration or widen C_FACTOR_RANGE.')