from pathlib import Pathimport csvimport jsonimport numpy as npimport matplotlibmatplotlib.use('Agg')import matplotlib.pyplot as pltimport matplotlib.tri as mtri
OUT = Path(__file__).resolve().parentscores = np.array([ [.88, .88, .70], [.90, .85, .71], [.92, .82, .72], [.91, .86, .69], [.89, .90, .67], [.87, .91, .68], [.86, .89, .71], [.70, .70, .60]])
def nondominated(a): return np.array([not np.any(np.all(a >= row, axis=1) & np.any(a > row, axis=1)) for row in a])
keep = nondominated(scores)nds = scores[keep]assert len(nds) == 7with (OUT/'synthetic_scores.csv').open('w', newline='', encoding='utf-8') as f: writer = csv.writer(f) writer.writerow(['solution_id', 'KGE_Q', 'KGE_Sstar', 'KGE_ET', 'is_NDS']) for k, row in enumerate(scores): writer.writerow([k+1, *row, int(keep[k])])
N = 100weights = np.array([[i, j, N-i-j] for i in range(N+1) for j in range(N+1-i)], float)/Nassert len(weights) == 5151assert np.all(weights >= 0) and np.allclose(weights.sum(1), 1)ce_each = nds @ weights.T ce = np.median(ce_each, axis=0)
def xy(w): w = np.asarray(w, float) return w[..., 1]+.5*w[..., 2], np.sqrt(3)/2*w[..., 2]
x, y = xy(weights)mesh = mtri.Triangulation(x, y)np.savetxt(OUT/'weight_surface.csv', np.column_stack([weights, x, y, ce]), delimiter=',', header='w_Q,w_Sstar,w_ET,x,y,CE_median', comments='')fig, ax = plt.subplots(figsize=(7.5, 6.4))levels = np.linspace(.60, .90, 13)filled = ax.tricontourf(mesh, ce, levels=levels, cmap='RdBu', extend='both')for level, color, style in [(.48, '#bdbdbd', '--'), (.65, '#bdbdbd', '-'), (.83, '#303030', '--')]: if ce.min() < level < ce.max(): lines = ax.tricontour(mesh, ce, levels=[level], colors=[color],
linestyles=[style], linewidths=1.8) ax.clabel(lines, fmt='%.2f', fontsize=10)for k in range(3): for t in [.2, .4, .6, .8]: p = np.zeros((2, 3)); p[:, k] = t p[0, (k+1)%3] = 1-t; p[1, (k+2)%3] = 1-t ax.plot(*xy(p), color='gray', lw=.55, ls='--', alpha=.65)h = np.sqrt(3)/2ax.plot([0, 1, .5, 0], [0, 0, h, 0], color='#333333', lw=1.4)for t in np.linspace(0, 1, 6): ax.text(t, -.035, f'{t:.1f}', ha='center', va='top', fontsize=10) ax.text(.5*(1-t)-.035, h*(1-t), f'{t:.1f}', ha='right', fontsize=10) ax.text(1-.5*t+.035, h*t, f'{t:.1f}', ha='left', fontsize=10)ax.text(.5, -.105, r'$w_{S^*}$', ha='center', fontsize=15)ax.text(.13, .45, r'$w_Q$', rotation=60, ha='center', fontsize=15)ax.text(.87, .45, r'$w_{ET}$', rotation=-60, ha='center', fontsize=15)ax.set(xlim=(-.14, 1.14), ylim=(-.13, h+.08), aspect='equal')ax.axis('off')ax.set_title('Panel b method reproduction | SYNTHETIC DATA', fontsize=13, pad=16)bar = fig.colorbar(filled, ax=ax, shrink=.8, pad=.06, ticks=[.60, .65, .70, .75, .80, .85, .90])bar.set_label('Median combined efficiency')fig.tight_layout()for ext in ['png', 'svg']: fig.savefig(OUT/f'panel_b_synthetic.{ext}', dpi=300, bbox_inches='tight')plt.close(fig)check_w = np.array([[1,0,0], [0,1,0], [0,0,1], [1/3]*3, [.6,.3,.1]])values = np.median(nds @ check_w.T, axis=0)assert np.allclose(values, [.89, .88, .70, .82, .866])report = dict(synthetic=True, candidate_count=len(scores), NDS_count=len(nds), weight_count=len(weights), grid_step=1/N, minimum=float(ce.min()), maximum=float(ce.max()), check_values=values.tolist())(OUT/'validation.json').write_text(json.dumps(report, indent=2), encoding=
'utf-8')print(json.dumps(report, indent=2))
