"""Read-only coverage and amplitude analysis of saved, phase-paired MQL8 curves. No simulation, re-pairing, resampling, filtering, or production changes. Integral metrics are trapezoidal estimates on the saved common grid, not bounds on any unsampled transient. Occupancy durations use sample-cell weights and are not located threshold-crossing times. Run with --plots in a matplotlib environment. """ from __future__ import annotations import argparse import hashlib import json from pathlib import Path import numpy as np ROOT = Path(__file__).resolve().parents[2] BASE = ROOT / 'test/lstp-mainline-20260917' def digest(path): return hashlib.sha256(path.read_bytes()).hexdigest() def ratio(numerator, denominator): return float(100 * numerator / denominator) if denominator else None def intervals(mask, time, weights): indices = np.flatnonzero(mask) blocks = np.split(indices, np.flatnonzero(np.diff(indices) > 1) + 1) return [dict(firstSample=float(time[b[0]]), lastSample=float(time[b[-1]]), sampleCount=len(b), cellDurationEstimate=float(weights[b].sum())) for b in blocks if len(b)] def mask_summary(mask, time, weights): return dict(sampleCount=int(mask.sum()), samplePercent=ratio(mask.sum(), len(mask)), cellDurationEstimate=float(weights[mask].sum()), timePercentEstimate=ratio(weights[mask].sum(), weights.sum()), intervals=intervals(mask, time, weights)) def curve_stats(row, actual, reference, time, weights): error = actual - reference absolute = np.abs(error) magnitude = np.abs(reference) epsilon = row['epsilon'] active = magnitude > epsilon near = ~active relative = np.zeros(len(time)) relative[active] = absolute[active] / magnitude[active] * 100 peak = float(magnitude.max()) significant = magnitude > max(epsilon, .01 * peak) near_bad = near & (absolute > epsilon) signed_integral = np.concatenate(([0.], np.cumsum( .5 * (error[1:] + error[:-1]) * np.diff(time)))) result = dict(key=row['key'], quantity=row['quantity'], unit=row['unit'], epsilon=epsilon, sampleCount=len(time), referencePeak=peak, activeCount=int(active.sum()), significantCount=int(significant.sum()), maximumAbsolute=float(absolute.max()), worstAbsoluteTime=float(time[np.argmax(absolute)]), maximumAbsolutePercentOfPeak=ratio(absolute.max(), peak), rmse=float(np.sqrt(np.dot(weights, error**2) / weights.sum())), relativeL2Percent=ratio(np.sqrt(np.dot(weights, error**2)), np.sqrt(np.dot(weights, reference**2))), integratedAbsoluteError=float(np.dot(weights, absolute)), integratedReferenceMagnitude=float(np.dot(weights, magnitude)), relativeL1Percent=ratio(np.dot(weights, absolute), np.dot(weights, magnitude)), signedIntegralError=float(signed_integral[-1]), maxCumulativeSignedError=float(np.abs(signed_integral).max()), activeRelativePercentiles={str(p): float(np.percentile(relative[active], p)) if active.any() else None for p in (50, 95, 99, 100)}, significantMaxRelativePercent=float(relative[significant].max()) if significant.any() else None, nearZeroCount=int(near.sum()), nearZeroMaxAbsolute=float(absolute[near].max()) if near.any() else None, nearZeroAboveEpsilon=mask_summary(near_bad, time, weights), sensitivity={str(factor): int(((magnitude > epsilon * factor) & (absolute > .05 * magnitude)).sum()) for factor in (.1, 1., 10.)}) masks = {} for threshold in (1, 5): mask = active & (relative > threshold) masks[str(threshold)] = mask result['above' + str(threshold)] = mask_summary(mask, time, weights) | dict( activePercent=ratio(mask.sum(), active.sum())) result['relative5OrNearZeroAbsolute'] = mask_summary(masks['5'] | near_bad, time, weights) result['examplesAbove5'] = [dict(time=float(time[i]), platform=float(actual[i]), amesim=float(reference[i]), absoluteError=float(absolute[i]), relativePercent=float(relative[i])) for i in np.flatnonzero(masks['5'])] result['nearZeroExamples'] = [dict(time=float(time[i]), platform=float(actual[i]), amesim=float(reference[i]), absoluteError=float(absolute[i])) for i in np.flatnonzero(near_bad)] # Cross-check the earlier diagnostic without changing its epsilon or pairing. assert result['above5']['sampleCount'] == row['above5PercentCount'] assert int(near_bad.sum()) == row['nearZeroBeyondEpsilon'] return result, masks['5'], near_bad def analyze(source): paths = [source / name for name in ('comparison.json', 'curves.npz')] before = {str(path): digest(path) for path in paths} comparison = json.loads(paths[0].read_bytes()) arrays = np.load(paths[1], allow_pickle=False) time = arrays['time'] assert arrays['phaseMatched'].all() and np.all(np.diff(time) > 0) dt = np.diff(time) weights = np.r_[dt[0] / 2, (dt[:-1] + dt[1:]) / 2, dt[-1] / 2] assert np.isclose(weights.sum(), time[-1] - time[0]) rows, masks, near_masks = [], {}, {} for row in comparison['curves']: actual, reference = (arrays[s + '|' + row['key']] for s in ('platform', 'amesim')) assert np.isfinite(actual).all() and np.isfinite(reference).all() stats, mask, near = curve_stats(row, actual, reference, time, weights) rows.append(stats) masks[row['key']], near_masks[row['key']] = mask, near groups = {} for quantity in sorted({r['quantity'] for r in rows}): selected = [r for r in rows if r['quantity'] == quantity] all_count = len(selected) * len(time) active_count = sum(r['activeCount'] for r in selected) union = np.any([masks[r['key']] for r in selected], axis=0) near_union = np.any([near_masks[r['key']] for r in selected], axis=0) def worst(field): candidates = [r for r in selected if r[field] is not None] if not candidates: return None r = max(candidates, key=lambda r: r[field]) return dict(key=r['key'], value=r[field]) groups[quantity] = dict(curveCount=len(selected), sampleCount=all_count, activeCount=active_count, above5Count=sum(r['above5']['sampleCount'] for r in selected), above5SamplePercent=ratio(sum(r['above5']['sampleCount'] for r in selected), all_count), above5ActivePercent=ratio(sum(r['above5']['sampleCount'] for r in selected), active_count), anyCurveAbove5=mask_summary(union, time, weights), nearZeroAboveEpsilonCount=sum(r['nearZeroAboveEpsilon']['sampleCount'] for r in selected), nearZeroAboveEpsilonSamplePercent=ratio(sum(r['nearZeroAboveEpsilon']['sampleCount'] for r in selected), all_count), anyCurveNearZeroAboveEpsilon=mask_summary(near_union, time, weights), anyCurveRelative5OrNearZeroAbsolute=mask_summary(union | near_union, time, weights), worst={field: worst(field) for field in ('maximumAbsolute', 'maximumAbsolutePercentOfPeak', 'relativeL2Percent', 'relativeL1Percent', 'significantMaxRelativePercent', 'integratedAbsoluteError', 'maxCumulativeSignedError')}) union = np.any(list(masks.values()), axis=0) near_union = np.any(list(near_masks.values()), axis=0) windows = [] for start, end in [(0., 1.), (1., 10.8), (10.8, 21.6), (21.6, 32.4), (32.4, 43.2), (43.2, 50.)]: include = (time >= start - 1e-12) & (time < end - 1e-12 if end < 50 else time <= end) details = {} for quantity in ('mass_flow', 'enthalpy_flow'): selected = [r for r in rows if r['quantity'] == quantity] details[quantity] = dict( above5Count=sum(int((masks[r['key']] & include).sum()) for r in selected), nearZeroAboveEpsilonCount=sum(int((near_masks[r['key']] & include).sum()) for r in selected), maximumAbsolute=max(float(np.abs(arrays['platform|' + r['key']] - arrays['amesim|' + r['key']])[include].max()) for r in selected)) windows.append(dict(start=start, end=end, sampleCount=int(include.sum()), groups=details)) assert sum(r['above5']['sampleCount'] for r in rows) == comparison['above5PercentCount'] assert all(digest(path) == before[str(path)] for path in paths) return dict(source=str(source), inputHashes=before, sourceUnchanged=True, grid=dict(start=float(time[0]), end=float(time[-1]), count=len(time), interval=float(np.median(dt))), curveCount=len(rows), affectedCurveCount=sum(r['above5']['sampleCount'] > 0 for r in rows), above5Count=sum(r['above5']['sampleCount'] for r in rows), anyCurveAbove5=mask_summary(union, time, weights), anyCurveNearZeroAboveEpsilon=mask_summary(near_union, time, weights), anyCurveRelative5OrNearZeroAbsolute=mask_summary(union | near_union, time, weights), groups=groups, windows=windows, curves=rows) def plots(source, out): import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib import font_manager font = Path('C:/Windows/Fonts/msyh.ttc') if font.exists(): font_manager.fontManager.addfont(str(font)) plt.rcParams['font.family'] = font_manager.FontProperties(fname=str(font)).get_name() plt.rcParams.update({'font.size': 10, 'axes.unicode_minus': False, 'axes.spines.top': False, 'axes.spines.right': False}) arrays = np.load(source / 'curves.npz', allow_pickle=False) t = arrays['time'] fig, axes = plt.subplots(2, 3, figsize=(15, 8), constrained_layout=True) selected = [('amesim_pn3node2_3.reference_mass_flow', '质量流量', 'kg/s', 1e6, 'mg/s'), ('amesim_p4node2_4.reference_enthalpy_flow', '焓流', 'W', 1., 'W')] for row, (key, label, unit, scale, small_unit) in enumerate(selected): y, ref = arrays['platform|' + key], arrays['amesim|' + key] for col, bounds in enumerate(((0, 50), (.27, .35))): ax = axes[row, col] mask = (t >= bounds[0] - 1e-12) & (t <= bounds[1] + 1e-12) factor = 1. if col == 0 else scale ax.plot(t[mask], ref[mask] * factor, color='#dd863b', lw=2, label='Amesim') ax.plot(t[mask], y[mask] * factor, color='#126ca6', lw=1, ls='--', label='平台') ax.set(xlabel='时间 / s', ylabel=unit if col == 0 else small_unit, title=label + ('：50 s 全程' if col == 0 else '：接近零的衰减尾部放大')) if col == 1: ax.axvspan(.295, .325, color='#d84b43', alpha=.12, label='差异集中区') ax.legend(fontsize=8) ax = axes[row, 2] group_keys = [k for k in arrays.files if k.startswith('platform|') and k.endswith('.reference_' + ('mass_flow' if row == 0 else 'enthalpy_flow'))] envelope = np.max([np.abs(arrays[k] - arrays[k.replace('platform|', 'amesim|', 1)]) for k in group_keys], axis=0) ax.plot(t, envelope * scale, color='#8b3d50', lw=.9) ax.set(xlabel='时间 / s', ylabel=small_unit, title=label + '：16 条曲线最大绝对差包络') for ax in axes[row]: ax.grid(alpha=.2) fig.suptitle('八路基线剩余差异：全程、初始衰减段与绝对差\n0–50 s，共同网格 10 ms；事件侧已配对；不代表网格间瞬态的误差上界', fontsize=14) fig.savefig(out / 'coverage.png', dpi=150) plt.close(fig) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--output', type=Path, default=ROOT / 'test/mql8-curve-coverage-20260917') parser.add_argument('--plots', action='store_true') args = parser.parse_args() args.output.mkdir(parents=True, exist_ok=True) sources = dict(cyclic=BASE / 'event-output-comparison', noncyclic=BASE / 'baseline/noncyclic') result = dict(method='Saved-grid statistics; original phase pairing and epsilon preserved. ' 'L1/L2 normalized by each reference curve; no time interpolation. ' 'Integrals and durations are grid estimates only. Noncyclic reference uses ordinary output.', profiles={name: analyze(path) for name, path in sources.items()}) (args.output / 'coverage.json').write_text(json.dumps(result, ensure_ascii=False, indent=2, allow_nan=False) + '\n', encoding='utf-8') if args.plots: plots(sources['cyclic'], args.output) for name, profile in result['profiles'].items(): print(name, 'above5:', profile['above5Count'], 'any time:', profile['anyCurveAbove5']) for quantity in ('force', 'mass_flow', 'enthalpy_flow', 'pressure', 'temperature', 'gap', 'velocity'): print(quantity, json.dumps(profile['groups'][quantity], ensure_ascii=False)) if __name__ == '__main__': main()