feat(SR): publish canonical analysis package
Canonicalize V5 case identities and preserve the SR evidence chain while replacing ambiguous diagnostics with reproducible tables, phase-matched flow fields, and presentation-ready summaries. Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
Frank14f
+433
@@ -0,0 +1,433 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Export a phase-matched Target/PPO/SR Kármán vorticity comparison.
|
||||
|
||||
CUDA-backed modules are imported only after argument validation. The exporter
|
||||
uses the canonical Stage-3 environment and policy constructors for controlled
|
||||
runs and mirrors ``build_karman_cloak_env`` exactly for the target-only run.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import hashlib
|
||||
import json
|
||||
import os
|
||||
import platform
|
||||
import subprocess
|
||||
import sys
|
||||
from datetime import datetime, timezone
|
||||
from pathlib import Path
|
||||
from typing import Any, Mapping, Sequence
|
||||
|
||||
import numpy as np
|
||||
|
||||
REPO_ROOT = Path(__file__).resolve().parents[3]
|
||||
SRC_ROOT = REPO_ROOT / "src"
|
||||
for root in (REPO_ROOT, SRC_ROOT):
|
||||
if str(root) not in sys.path:
|
||||
sys.path.insert(0, str(root))
|
||||
|
||||
from SR_analysis.configs import FIFO_LEN, LEGACY_CFG_DIR, get_scene # noqa: E402
|
||||
from SR_analysis.stage_3_validate import ( # noqa: E402
|
||||
DATA_TYPE,
|
||||
ValidationPlan,
|
||||
build_karman_environment,
|
||||
build_policy,
|
||||
load_formula_pair,
|
||||
prepare_plan,
|
||||
)
|
||||
from SR_analysis.utils.provenance import atomic_write_json, hash_file, hash_json # noqa: E402
|
||||
|
||||
STEM = "07_flow_field_comparison_karman_re100"
|
||||
SCHEMA = "sr-flow-comparison-v2"
|
||||
CANDIDATE_DTYPE = np.dtype(np.float32)
|
||||
DEFAULT_PACKAGE = REPO_ROOT / "src/SR_analysis/results/runs/article2-plotting-package-20260721"
|
||||
DEFAULT_FORMULAS = REPO_ROOT / "src/SR_analysis/results/runs/article-refit-karman-topology-a-20260718/formulas"
|
||||
DEFAULT_ALIGNMENT = DEFAULT_PACKAGE / "phase_alignment.json"
|
||||
D_LATTICE = 20.0
|
||||
SENSOR_LAYOUT = ("upper_ux", "upper_uy", "center_ux", "center_uy", "lower_ux", "lower_uy")
|
||||
|
||||
|
||||
def sha256_bytes(data: bytes) -> str:
|
||||
return hashlib.sha256(data).hexdigest()
|
||||
|
||||
|
||||
def git_sha() -> str | None:
|
||||
try:
|
||||
return subprocess.run(["git", "-C", str(REPO_ROOT), "rev-parse", "HEAD"], check=True, capture_output=True, text=True).stdout.strip()
|
||||
except (OSError, subprocess.CalledProcessError):
|
||||
return None
|
||||
|
||||
|
||||
def wrapped_phase_difference(a: np.ndarray | float, b: np.ndarray | float) -> np.ndarray:
|
||||
"""Signed shortest angular difference ``a-b`` in [-pi, pi]."""
|
||||
delta = np.asarray(a, dtype=np.float64) - np.asarray(b, dtype=np.float64)
|
||||
return np.arctan2(np.sin(delta), np.cos(delta))
|
||||
|
||||
|
||||
def standardized_center_phase(trace: np.ndarray, start: int, stop: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
|
||||
"""Return center-sensor phase from per-trajectory stable-window z scores."""
|
||||
sensors = np.asarray(trace, dtype=np.float64)
|
||||
if sensors.ndim != 2 or sensors.shape[1] < 4 or not 0 <= start < stop <= len(sensors):
|
||||
raise ValueError("trace/window must provide center (ux,uy) channels")
|
||||
center = sensors[start:stop, 2:4]
|
||||
mean = center.mean(axis=0)
|
||||
scale = center.std(axis=0, ddof=0)
|
||||
if np.any(~np.isfinite(scale)) or np.any(scale <= np.finfo(float).eps):
|
||||
raise ValueError("stable-window center-sensor scale is zero or non-finite")
|
||||
z = (center - mean) / scale
|
||||
return np.arctan2(z[:, 1], z[:, 0]), mean, scale
|
||||
|
||||
|
||||
def select_joint_phase_match(
|
||||
target_trace: np.ndarray,
|
||||
ppo_trace: np.ndarray,
|
||||
sr_trace: np.ndarray,
|
||||
start: int,
|
||||
stop: int,
|
||||
*,
|
||||
temporal_weight: float = 0.001,
|
||||
) -> dict[str, Any]:
|
||||
"""Jointly match downstream center-sensor limit-cycle phase.
|
||||
|
||||
Exhaustive stable-window search minimizes the two wrapped phase errors plus
|
||||
``temporal_weight * (|i_ppo-i_target| + |i_sr-i_target|)``. Requiring the
|
||||
same local phase direction rejects branch-reversed phase-portrait matches.
|
||||
"""
|
||||
if temporal_weight < 0:
|
||||
raise ValueError("temporal_weight must be non-negative")
|
||||
phases, means, scales, directions = [], {}, {}, []
|
||||
for label, trace in (("Target", target_trace), ("PPO", ppo_trace), ("SR", sr_trace)):
|
||||
theta, mean, scale = standardized_center_phase(trace, start, stop)
|
||||
phases.append(theta); means[label] = mean.tolist(); scales[label] = scale.tolist()
|
||||
directions.append(np.sign(np.gradient(np.unwrap(theta))))
|
||||
best = None
|
||||
for ti, target_index in enumerate(range(start, stop)):
|
||||
for pi, ppo_index in enumerate(range(start, stop)):
|
||||
if directions[1][pi] != directions[0][ti]:
|
||||
continue
|
||||
for si, sr_index in enumerate(range(start, stop)):
|
||||
if directions[2][si] != directions[0][ti]:
|
||||
continue
|
||||
ppo_error = abs(float(wrapped_phase_difference(phases[1][pi], phases[0][ti])))
|
||||
sr_error = abs(float(wrapped_phase_difference(phases[2][si], phases[0][ti])))
|
||||
temporal_distance = abs(ppo_index-target_index) + abs(sr_index-target_index)
|
||||
penalty = temporal_weight * temporal_distance
|
||||
objective = ppo_error + sr_error + penalty
|
||||
key = (objective, temporal_distance, target_index, ppo_index, sr_index)
|
||||
if best is None or key < best[0]:
|
||||
best = (key, ti, pi, si, ppo_error, sr_error, penalty)
|
||||
if best is None:
|
||||
raise ValueError("no same-direction phase match in stable window")
|
||||
key, ti, pi, si, ppo_error, sr_error, penalty = best
|
||||
indices = {"Target": start+ti, "PPO": start+pi, "SR": start+si}
|
||||
return {
|
||||
"indices": indices,
|
||||
"phase_angles_rad": {"Target": float(phases[0][ti]), "PPO": float(phases[1][pi]), "SR": float(phases[2][si])},
|
||||
"wrapped_angle_errors_rad": {"PPO": ppo_error, "SR": sr_error},
|
||||
"local_phase_direction": {"Target": int(directions[0][ti]), "PPO": int(directions[1][pi]), "SR": int(directions[2][si])},
|
||||
"stable_window_center_mean": means,
|
||||
"stable_window_center_std": scales,
|
||||
"objective_terms": {"angle_error_sum_rad": ppo_error+sr_error, "temporal_distance_samples": int(key[1]), "temporal_weight_rad_per_sample": temporal_weight, "temporal_penalty_rad": penalty, "objective": float(key[0])},
|
||||
"candidate_count_per_trajectory": stop-start,
|
||||
}
|
||||
|
||||
|
||||
def orient_vorticity_xy_to_yx(omega_xy: np.ndarray, field_shape: Sequence[int]) -> np.ndarray:
|
||||
"""Convert vorticity_from_ddf's (NX, NY) result to image (NY, NX)."""
|
||||
nx, ny = map(int, field_shape[:2])
|
||||
omega = np.asarray(omega_xy)
|
||||
if omega.shape != (nx, ny):
|
||||
raise ValueError(f"expected vorticity shape {(nx, ny)}, got {omega.shape}")
|
||||
return omega.T.copy()
|
||||
|
||||
|
||||
def crop_yx(field_yx: np.ndarray, xlim_d: Sequence[float], ylim_d: Sequence[float], *, d_lattice: float = D_LATTICE, center_y_lattice: float | None = None) -> tuple[np.ndarray, dict[str, Any]]:
|
||||
"""Crop an image-oriented field using physical x/D and centered y/D."""
|
||||
field = np.asarray(field_yx)
|
||||
if field.ndim != 2:
|
||||
raise ValueError("field must be 2-D in (y, x) order")
|
||||
ny, nx = field.shape
|
||||
cy = (ny - 1) / 2 if center_y_lattice is None else float(center_y_lattice)
|
||||
x0 = max(0, int(np.ceil(float(xlim_d[0]) * d_lattice)))
|
||||
x1 = min(nx, int(np.floor(float(xlim_d[1]) * d_lattice)) + 1)
|
||||
y0 = max(0, int(np.ceil(cy + float(ylim_d[0]) * d_lattice)))
|
||||
y1 = min(ny, int(np.floor(cy + float(ylim_d[1]) * d_lattice)) + 1)
|
||||
if x0 >= x1 or y0 >= y1:
|
||||
raise ValueError("requested crop does not intersect field")
|
||||
extent = ((x0 / d_lattice), ((x1 - 1) / d_lattice), ((y0 - cy) / d_lattice), ((y1 - 1 - cy) / d_lattice))
|
||||
return field[y0:y1, x0:x1].copy(), {"x_slice": [x0, x1], "y_slice": [y0, y1], "extent_xD_yD": list(extent)}
|
||||
|
||||
|
||||
def manifest_artifacts(package_dir: Path, repo_root: Path = REPO_ROOT) -> list[dict[str, str]]:
|
||||
suffixes = {".csv", ".png", ".pdf", ".md", ".json", ".npz"}
|
||||
return [
|
||||
{"path": str(path.relative_to(repo_root)), "sha256": hash_file(path)}
|
||||
for path in sorted(package_dir.rglob("*"))
|
||||
if path.is_file() and path.suffix.lower() in suffixes and path.name != "manifest.json"
|
||||
]
|
||||
|
||||
|
||||
def update_package_manifest(package_dir: Path, repo_root: Path = REPO_ROOT) -> None:
|
||||
path = package_dir / "manifest.json"
|
||||
old = json.loads(path.read_text(encoding="utf-8")) if path.is_file() else {}
|
||||
summary = dict(old.get("summary", {}))
|
||||
summary["publication_figures"] = max(5, int(summary.get("publication_figures", 0)))
|
||||
summary["flow_field_npz"] = 1
|
||||
summary["presentation_pages"] = 2
|
||||
atomic_write_json(path, {"schema_version": "sr-plotting-package-v3", "source_policy": old.get("source_policy", "immutable article artifacts; no scientific refit"), "summary": summary, "artifacts": manifest_artifacts(package_dir, repo_root)})
|
||||
|
||||
|
||||
def _runtime_cfd() -> tuple[Any, Any, Any]:
|
||||
from LegacyCelerisLab import FlowField
|
||||
from SR_analysis.utils.cfd_interface import load_legacy_configs, vorticity_from_ddf
|
||||
return FlowField, load_legacy_configs, vorticity_from_ddf
|
||||
|
||||
|
||||
def _candidate_index_array(start: int, stop: int) -> np.ndarray:
|
||||
"""Return the exact trace indices represented by stable-window fields."""
|
||||
if not 0 <= start < stop:
|
||||
raise ValueError("candidate window must satisfy 0 <= start < stop")
|
||||
return np.arange(start, stop, dtype=np.int64)
|
||||
|
||||
|
||||
def select_exact_candidate_field(
|
||||
candidates: np.ndarray,
|
||||
candidate_trace_indices: np.ndarray,
|
||||
selected_trace_index: int,
|
||||
*,
|
||||
expected_start: int,
|
||||
expected_stop: int,
|
||||
) -> tuple[np.ndarray, int]:
|
||||
"""Select a field by its same-run trace index, rejecting mapping drift."""
|
||||
fields = np.asarray(candidates)
|
||||
indices = np.asarray(candidate_trace_indices)
|
||||
expected = _candidate_index_array(expected_start, expected_stop)
|
||||
if fields.ndim != 3 or fields.dtype != CANDIDATE_DTYPE:
|
||||
raise ValueError("candidate fields must be a float32 (sample,y,x) array")
|
||||
if indices.ndim != 1 or not np.issubdtype(indices.dtype, np.integer):
|
||||
raise ValueError("candidate trace indices must be a one-dimensional integer array")
|
||||
if fields.shape[0] != len(indices):
|
||||
raise ValueError("candidate field/index counts differ")
|
||||
if not np.array_equal(indices, expected):
|
||||
raise ValueError("candidate trace-index mapping does not exactly cover the stable window")
|
||||
matches = np.flatnonzero(indices == int(selected_trace_index))
|
||||
if len(matches) != 1:
|
||||
raise IndexError("selected trace index has no unique same-run candidate field")
|
||||
candidate_index = int(matches[0])
|
||||
if candidate_index != int(selected_trace_index) - expected_start:
|
||||
raise AssertionError("candidate offset and trace index disagree")
|
||||
return fields[candidate_index], candidate_index
|
||||
|
||||
|
||||
def _allocate_candidates(field_shape: Sequence[int], start: int, stop: int) -> np.ndarray:
|
||||
nx, ny = map(int, field_shape[:2])
|
||||
return np.empty((stop - start, ny, nx), dtype=CANDIDATE_DTYPE)
|
||||
|
||||
|
||||
def _store_candidate(
|
||||
candidates: np.ndarray,
|
||||
trace_index: int,
|
||||
start: int,
|
||||
omega_xy: np.ndarray,
|
||||
field_shape: Sequence[int],
|
||||
) -> None:
|
||||
candidate_index = trace_index - start
|
||||
if not 0 <= candidate_index < len(candidates):
|
||||
raise IndexError("candidate trace index lies outside allocated stable window")
|
||||
field = D_LATTICE * orient_vorticity_xy_to_yx(omega_xy, field_shape)
|
||||
candidates[candidate_index] = np.asarray(field, dtype=CANDIDATE_DTYPE)
|
||||
|
||||
|
||||
def _target_trace_and_candidates(
|
||||
cfg: Mapping[str, Any], device: int, n_samples: int, candidate_start: int
|
||||
) -> tuple[np.ndarray, np.ndarray, np.ndarray, tuple[int, int]]:
|
||||
"""Run Target once and couple every stable trace sample to its field."""
|
||||
FlowField, load_configs, vorticity_from_ddf = _runtime_cfd()
|
||||
cuda_cfg, field_cfg = load_configs(LEGACY_CFG_DIR)
|
||||
ff = FlowField(field_cfg._replace(viscosity=float(cfg["nu"])), cuda_cfg, device_id=device)
|
||||
try:
|
||||
cy = (ff.FIELD_SHAPE[1] - 1) / 2
|
||||
ff.add_cylinder((10.0 * D_LATTICE, cy, 0.0), D_LATTICE)
|
||||
for y_off in (2.0, 0.0, -2.0):
|
||||
ff.add_sensor((40.0 * D_LATTICE, cy + y_off * D_LATTICE, 0.0), D_LATTICE / 4.0)
|
||||
n_obj = ff.obs.size // 2
|
||||
zero = np.zeros(n_obj, dtype=DATA_TYPE)
|
||||
ff.run(int(4 * ff.FIELD_SHAPE[0] / float(cfg["u0"])), zero)
|
||||
rows: list[np.ndarray] = []
|
||||
candidates = _allocate_candidates(ff.FIELD_SHAPE, candidate_start, n_samples)
|
||||
candidate_indices = _candidate_index_array(candidate_start, n_samples)
|
||||
for index in range(n_samples):
|
||||
ff.run(int(cfg["sample_interval"]), zero)
|
||||
rows.append(ff.obs.copy()[2:8].astype(np.float64))
|
||||
if index >= candidate_start:
|
||||
_store_candidate(candidates, index, candidate_start, vorticity_from_ddf(ff, float(cfg["u0"])), ff.FIELD_SHAPE)
|
||||
if len(rows) != n_samples or len(candidates) != n_samples - candidate_start:
|
||||
raise AssertionError("Target trace/candidate collection is incomplete")
|
||||
return np.asarray(rows), candidates, candidate_indices, tuple(map(int, ff.FIELD_SHAPE[:2]))
|
||||
finally:
|
||||
del ff
|
||||
|
||||
|
||||
def _controlled_trace_and_candidates(
|
||||
plan: ValidationPlan, device: int, n_samples: int, candidate_start: int
|
||||
) -> tuple[np.ndarray, np.ndarray, np.ndarray, tuple[int, int]]:
|
||||
"""Run one controller once and couple stable trace samples to fields."""
|
||||
env = build_karman_environment(plan, device)
|
||||
try:
|
||||
policy = build_policy(plan)
|
||||
raw = np.asarray(env.current_raw, dtype=np.float64)
|
||||
rows: list[np.ndarray] = []
|
||||
candidates = _allocate_candidates(env.ff.FIELD_SHAPE, candidate_start, n_samples)
|
||||
candidate_indices = _candidate_index_array(candidate_start, n_samples)
|
||||
_, _, vorticity_from_ddf = _runtime_cfd()
|
||||
for index in range(n_samples):
|
||||
omega, _, _ = policy.action(raw, index)
|
||||
raw = env.step(omega)
|
||||
policy.observe(omega)
|
||||
rows.append(raw[:6].copy())
|
||||
if index >= candidate_start:
|
||||
_store_candidate(candidates, index, candidate_start, vorticity_from_ddf(env.ff, float(plan.cfg["u0"])), env.ff.FIELD_SHAPE)
|
||||
if len(rows) != n_samples or len(candidates) != n_samples - candidate_start:
|
||||
raise AssertionError("controlled trace/candidate collection is incomplete")
|
||||
return np.asarray(rows), candidates, candidate_indices, tuple(map(int, env.ff.FIELD_SHAPE[:2]))
|
||||
finally:
|
||||
env.close()
|
||||
|
||||
|
||||
def _make_plan(scene: str, mode: str, n_steps: int, pair: Any | None, model_device: str) -> ValidationPlan:
|
||||
plan = prepare_plan(scene=scene, mode=mode, n_steps=n_steps, run_id="flow-comparison-in-memory", output_root=REPO_ROOT / ".flow-comparison-unused", formula_pair=pair)
|
||||
return ValidationPlan(**{**plan.__dict__, "cfg": {**plan.cfg, "model_device": model_device}})
|
||||
|
||||
|
||||
def _plot(fields: Mapping[str, np.ndarray], crop_meta: Mapping[str, Any], cfg: Mapping[str, Any], output_dir: Path, vmax: float) -> None:
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib.patches import Circle
|
||||
|
||||
fig, axes = plt.subplots(1, 3, figsize=(14.5, 4.3), sharex=True, sharey=True, constrained_layout=True)
|
||||
extent = crop_meta["extent_xD_yD"]
|
||||
for ax, label in zip(axes, ("Target", "PPO", "SR")):
|
||||
ax.imshow(fields[label], origin="lower", extent=extent, aspect="equal", cmap="RdBu_r", vmin=-vmax, vmax=vmax, interpolation="nearest")
|
||||
geometry = [(10.0, 0.0, 1.0, "disturbance")] if label == "Target" else [
|
||||
(10.0, 0.0, 1.0, "disturbance"),
|
||||
(float(cfg["pinball_front_x"]), 0.0, 0.5, "front"),
|
||||
(float(cfg["pinball_rear_x"]), 0.75, 0.5, "upper"),
|
||||
(float(cfg["pinball_rear_x"]), -0.75, 0.5, "lower"),
|
||||
]
|
||||
for x, y, radius, name in geometry:
|
||||
ax.add_patch(Circle((x, y), radius, facecolor="white", edgecolor="black", linewidth=0.8, zorder=4))
|
||||
ax.scatter([40.0] * 3, [2.0, 0.0, -2.0], s=12, marker="x", color="black", linewidths=0.8, zorder=5)
|
||||
ax.set_title(label)
|
||||
ax.set_xlabel(r"$x/D$")
|
||||
axes[0].set_ylabel(r"$y/D$")
|
||||
for suffix in ("png", "pdf"):
|
||||
fig.savefig(output_dir / f"{STEM}.{suffix}", dpi=300 if suffix == "png" else None)
|
||||
plt.close(fig)
|
||||
|
||||
|
||||
def build_parser() -> argparse.ArgumentParser:
|
||||
parser = argparse.ArgumentParser(description=__doc__)
|
||||
parser.add_argument("--scene", default="karman_re100", choices=("karman_re100",), help="explicit canonical scene")
|
||||
parser.add_argument("--package-dir", type=Path, default=DEFAULT_PACKAGE)
|
||||
parser.add_argument("--alignment-metadata", type=Path, default=DEFAULT_ALIGNMENT)
|
||||
parser.add_argument("--formula-front", type=Path, default=DEFAULT_FORMULAS / "joint_front.json")
|
||||
parser.add_argument("--formula-rear", type=Path, default=DEFAULT_FORMULAS / "joint_rear_shared_upper.json")
|
||||
parser.add_argument("--temporal-weight", type=float, default=0.001, help="phase objective penalty in radians per sample of temporal separation")
|
||||
parser.add_argument("--device", type=int, default=0, help="logical CFD device after CUDA_VISIBLE_DEVICES masking")
|
||||
parser.add_argument("--model-device", choices=("cpu",), default="cpu")
|
||||
parser.add_argument("--xlim", nargs=2, type=float, default=(7.0, 48.0), metavar=("XMIN_D", "XMAX_D"))
|
||||
parser.add_argument("--ylim", nargs=2, type=float, default=(-6.0, 6.0), metavar=("YMIN_D", "YMAX_D"))
|
||||
parser.add_argument("--replace", action="store_true")
|
||||
return parser
|
||||
|
||||
|
||||
def main(argv: Sequence[str] | None = None) -> int:
|
||||
args = build_parser().parse_args(argv)
|
||||
package = args.package_dir.resolve()
|
||||
outputs = [package / f"{STEM}.npz", package / f"{STEM}.json", package / "figures" / f"{STEM}.png", package / "figures" / f"{STEM}.pdf"]
|
||||
existing = [p for p in outputs if p.exists()]
|
||||
if existing and not args.replace:
|
||||
raise FileExistsError("refusing to overwrite: " + ", ".join(map(str, existing)))
|
||||
alignment = json.loads(args.alignment_metadata.read_text(encoding="utf-8"))
|
||||
start, stop = int(alignment["target_index_start_inclusive"]), int(alignment["target_index_stop_exclusive"])
|
||||
cfg = get_scene(args.scene)
|
||||
if cfg["scene_id"] != "karman":
|
||||
raise ValueError("only the Kármán environment is supported")
|
||||
pair = load_formula_pair(args.formula_front, args.formula_rear)
|
||||
n_samples = stop
|
||||
ppo_plan = _make_plan(args.scene, "ppo", n_samples, None, args.model_device)
|
||||
sr_plan = _make_plan(args.scene, "pysr", n_samples, pair, args.model_device)
|
||||
|
||||
# Each condition is initialized and run exactly once. The three float32
|
||||
# stable-window buffers remain in host memory until joint phase selection.
|
||||
target_trace, target_candidates, target_candidate_indices, shape = _target_trace_and_candidates(
|
||||
cfg, args.device, n_samples, start
|
||||
)
|
||||
ppo_trace, ppo_candidates, ppo_candidate_indices, ppo_shape = _controlled_trace_and_candidates(
|
||||
ppo_plan, args.device, n_samples, start
|
||||
)
|
||||
sr_trace, sr_candidates, sr_candidate_indices, sr_shape = _controlled_trace_and_candidates(
|
||||
sr_plan, args.device, n_samples, start
|
||||
)
|
||||
if shape != ppo_shape or shape != sr_shape:
|
||||
raise ValueError("field shapes differ between cases")
|
||||
|
||||
phase_match = select_joint_phase_match(target_trace, ppo_trace, sr_trace, start, stop, temporal_weight=args.temporal_weight)
|
||||
selected_indices = phase_match["indices"]
|
||||
candidate_sets = {
|
||||
"Target": (target_candidates, target_candidate_indices),
|
||||
"PPO": (ppo_candidates, ppo_candidate_indices),
|
||||
"SR": (sr_candidates, sr_candidate_indices),
|
||||
}
|
||||
full: dict[str, np.ndarray] = {}
|
||||
selected_candidate_indices: dict[str, int] = {}
|
||||
for label, (candidates, trace_indices) in candidate_sets.items():
|
||||
field, candidate_index = select_exact_candidate_field(
|
||||
candidates, trace_indices, selected_indices[label], expected_start=start, expected_stop=stop
|
||||
)
|
||||
full[label] = field.copy()
|
||||
selected_candidate_indices[label] = candidate_index
|
||||
candidate_field_bytes = int(sum(candidates.nbytes for candidates, _ in candidate_sets.values()))
|
||||
cropped: dict[str, np.ndarray] = {}
|
||||
crop_meta = None
|
||||
for label, field in full.items():
|
||||
cropped[label], meta = crop_yx(field, args.xlim, args.ylim)
|
||||
crop_meta = crop_meta or meta
|
||||
if meta != crop_meta:
|
||||
raise AssertionError("inconsistent crop metadata")
|
||||
finite = np.concatenate([np.abs(value[np.isfinite(value)]) for value in cropped.values()])
|
||||
if finite.size == 0:
|
||||
raise FloatingPointError("captured fields contain no finite vorticity")
|
||||
vmax = float(np.percentile(finite, 99.5))
|
||||
if not np.isfinite(vmax) or vmax <= 0:
|
||||
raise FloatingPointError("invalid shared color normalization")
|
||||
|
||||
package.mkdir(parents=True, exist_ok=True)
|
||||
(package / "figures").mkdir(parents=True, exist_ok=True)
|
||||
npz_path = outputs[0]
|
||||
np.savez_compressed(npz_path, target_vorticity_yx=full["Target"], ppo_vorticity_yx=full["PPO"], sr_vorticity_yx=full["SR"], target_sensors=target_trace, ppo_sensors=ppo_trace, sr_sensors=sr_trace, selected_indices=np.asarray([selected_indices[x] for x in ("Target", "PPO", "SR")]), phase_angles_rad=np.asarray([phase_match["phase_angles_rad"][x] for x in ("Target", "PPO", "SR")]), wrapped_angle_errors_rad=np.asarray([0.0, phase_match["wrapped_angle_errors_rad"]["PPO"], phase_match["wrapped_angle_errors_rad"]["SR"]]), local_phase_direction=np.asarray([phase_match["local_phase_direction"][x] for x in ("Target", "PPO", "SR")]), selected_candidate_indices=np.asarray([selected_candidate_indices[x] for x in ("Target", "PPO", "SR")]), extent_xD_yD=np.asarray(crop_meta["extent_xD_yD"]), crop_x_slice=np.asarray(crop_meta["x_slice"]), crop_y_slice=np.asarray(crop_meta["y_slice"]), field_shape_xy=np.asarray(shape), sensor_layout=np.asarray(SENSOR_LAYOUT))
|
||||
_plot(cropped, crop_meta, cfg, package / "figures", vmax)
|
||||
|
||||
dt = float(cfg["control_dt"])
|
||||
model_path = Path(ppo_plan.model_path) if ppo_plan.model_path else None
|
||||
metadata = {
|
||||
"schema_version": SCHEMA, "scene": args.scene, "description": "single deterministic downstream center-sensor limit-cycle phase snapshot; not full-field identity or an ensemble",
|
||||
"phase_matching": {"method": "one run per condition with same-sample stable-window field capture; exhaustive joint phase search", "phase_definition": "theta=atan2(z(center_uy), z(center_ux)); each trajectory standardized separately over stable window", "objective_formula": "|wrap(theta_PPO-theta_Target)| + |wrap(theta_SR-theta_Target)| + w*(|i_PPO-i_Target|+|i_SR-i_Target|)", "target_window": [start, stop], "direction_check": "all selected local unwrapped-phase derivatives have the same sign", "figure06_alignment_context_only": alignment, **phase_match},
|
||||
"selected": {label: {"index": int(selected_indices[label]), "candidate_index": int(selected_candidate_indices[label]), "t_D_over_U0": float(selected_indices[label] * dt), "phase_angle_rad": float(phase_match["phase_angles_rad"][label]), "wrapped_phase_error_rad": 0.0 if label == "Target" else float(phase_match["wrapped_angle_errors_rad"][label])} for label in ("Target", "PPO", "SR")},
|
||||
"sample_coupling": {"contract": "for every condition, the selected field and phase-diagnostic sensor values come from the same CFD sample in the same run", "candidate_trace_index_mapping": "candidate_index = trace_index - stable_window_start; exact contiguous mapping asserted before selection", "candidate_window": [start, stop], "candidate_count_per_condition": stop-start, "candidate_dtype": CANDIDATE_DTYPE.name, "candidate_field_shape_yx": [int(shape[1]), int(shape[0])], "candidate_host_memory_bytes": candidate_field_bytes, "candidate_host_memory_mib": candidate_field_bytes / (1024**2), "disk_contract": "only the three selected full fields and complete sensor traces are stored; stable-window candidates are memory-only"},
|
||||
"field": {"quantity": "omega_z D/U0", "source": "vorticity_from_ddf", "source_shape_order": "(NX,NY)", "stored_shape_order": "(NY,NX)", "crop": crop_meta, "shared_symmetric_vmax_percentile": {"percentile": 99.5, "vmax": vmax}, "geometry_D": {"disturbance": [10.0, 0.0, 1.0], "pinball": [[cfg["pinball_front_x"], 0.0, 0.5], [cfg["pinball_rear_x"], 0.75, 0.5], [cfg["pinball_rear_x"], -0.75, 0.5]], "sensors": [[40.0, 2.0], [40.0, 0.0], [40.0, -2.0]]}},
|
||||
"contracts": {"target": "exact build_karman_cloak_env geometry, stabilization, and sample interval; target has no pinball", "controlled": "stage_3_validate.prepare_plan/build_karman_environment/build_policy", "ppo_model_device": args.model_device, "cfd_logical_device": args.device},
|
||||
"hashes": {"config": hash_json(cfg), "formula_front": hash_file(args.formula_front), "formula_rear": hash_file(args.formula_rear), "formula_pair": pair.pair_hash, "ppo_model": hash_file(model_path) if model_path and model_path.is_file() else None, "alignment_metadata": hash_file(args.alignment_metadata), "npz": hash_file(npz_path), "png": hash_file(outputs[2]), "pdf": hash_file(outputs[3]), "exporter": hash_file(Path(__file__))},
|
||||
"paths": {"npz": str(npz_path.relative_to(REPO_ROOT)), "formula_front": str(args.formula_front.resolve()), "formula_rear": str(args.formula_rear.resolve()), "ppo_model": str(model_path) if model_path else None},
|
||||
"provenance": {"git_sha": git_sha(), "command": " ".join(sys.argv), "created_utc": datetime.now(timezone.utc).isoformat(), "python": sys.version, "platform": platform.platform(), "numpy": np.__version__, "CUDA_VISIBLE_DEVICES": os.environ.get("CUDA_VISIBLE_DEVICES")},
|
||||
}
|
||||
metadata["record_hash"] = hash_json(metadata)
|
||||
atomic_write_json(outputs[1], metadata)
|
||||
update_package_manifest(package)
|
||||
print(json.dumps({"outputs": [str(p) for p in outputs], "selected": metadata["selected"], "vmax": vmax}, indent=2))
|
||||
return 0
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main())
|
||||
Reference in New Issue
Block a user