{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"abe2d26a","cell_type":"markdown","source":"\n# TBD — uniform layout selection and PSNR-interval analysis\n\nThis notebook is the submission run for **SZ3** and **HPEZ** on five datasets:\nAMS/GEFS, longitudinal MR-MS, INSTANT-ODC, HyperLeaf2024, and IXI.\n\n**Selection rule.** The same procedure is used for every dataset–compressor pair:\n\n1. enumerate every semantic-axis permutation and contiguous folding supported by the 4-D codec interface;\n2. compute the locality statistic over the complete layout space;\n3. retain a uniform shortlist consisting of the strongest global layouts and the strongest layouts within **every folding partition**;\n4. evaluate that shortlist on development data;\n5. select one layout from disjoint validation data using measured PSNR and bit rate;\n6. evaluate only the selected and reference layouts on held-out statistical units.\n\nThe native-order reference layout is retained only as the comparison baseline.\n\n**Measured-PSNR intervals.** After all test curves are complete, the notebook applies the same six 10-dB measured-PSNR intervals to every dataset and compressor. This analysis is descriptive and never participates in layout selection. It reports per-unit values, dataset means with bootstrap confidence intervals, and an automatic interval-shape classification.\n","metadata":{}},{"id":"67114a8e","cell_type":"code","source":"\n# Execution controls only. These do not specify or favor any layout.\nimport os\nos.environ[\"HOC_FRESH\"] = \"1\" # use \"0\" only to resume an interrupted compatible run\nos.environ[\"HOC_RESUME\"] = \"1\"\nos.environ[\"HOC_STRICT_DATASETS\"] = \"1\"\nos.environ[\"HOC_STRICT_CODECS\"] = \"1\"\nos.environ[\"HOC_STRICT_ALL_CODEC_GATE\"] = \"1\"\nos.environ[\"HOC_WORK\"] = \"/kaggle/working/tbd_dimension_layout_final_v2\"\nos.environ[\"HOC_INPUT_ROOT\"] = \"/kaggle/input\"\n\n# Uniform search-coverage constants. They apply to every dataset and both codecs.\nos.environ[\"HOC_PROXY_GLOBAL_KEEP\"] = \"24\"\nos.environ[\"HOC_PROXY_PER_PARTITION_KEEP\"] = \"48\"\nos.environ[\"HOC_COARSE_KEEP\"] = \"24\"\n\n# Uniform post-test measured-PSNR reporting grid. It is not used for selection.\nos.environ[\"HOC_PSNR_INTERVAL_START_DB\"] = \"65\"\nos.environ[\"HOC_PSNR_INTERVAL_WIDTH_DB\"] = \"10\"\nos.environ[\"HOC_PSNR_INTERVAL_COUNT\"] = \"6\"\nos.environ[\"HOC_PSNR_INTERVAL_MIN_COVERAGE_DB\"] = \"5\"\n\nprint(\"Final TBD protocol configured: no dataset-specific layout identities; five datasets; SZ3 + HPEZ.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-19T00:40:43.82222Z","iopub.execute_input":"2026-09-19T00:40:43.823064Z","iopub.status.idle":"2026-09-19T00:40:43.839485Z","shell.execute_reply.started":"2026-09-19T00:40:43.823005Z","shell.execute_reply":"2026-09-19T00:40:43.838472Z"}},"outputs":[],"execution_count":null},{"id":"5978cb7b","cell_type":"markdown","source":"## Execute the benchmark","metadata":{}},{"id":"bb49ddc3","cell_type":"code","source":"#!/usr/bin/env python3\n\"\"\"Final five-domain benchmark for validation-selected dimension permutation and folding.\n\nDatasets evaluated under one protocol:\n * AMS/GEFS: variable x date x ensemble x lead x y x x (6-D)\n * Longitudinal MR-MS: patient x study x modality x z x y x x (6-D)\n * INSTANT-ODC: subject x modality x z x y x x (5-D)\n * HyperLeaf2024: plot x leaf x wavelength x y x x (5-D)\n * Preprocessed IXI MRI: subject x map x z x y x x (5-D)\n\nCompressors:\n * SZ3 (CPU, 4-D interface)\n * HPEZ/QoZ2 (CPU, 4-D interface)\n\nThe protocol enumerates semantic-axis permutations and contiguous mixed-radix foldings,\nuses a uniform global-plus-fold-partition shortlist on development data, selects one\nrepresentation per dataset and codec on disjoint validation data, and evaluates that\nrepresentation on independent held-out units. No dataset-specific layout identity is injected. Reference and selected\nrepresentations use the same absolute-tolerance grid. Every retained reconstruction\nmust satisfy the common pointwise absolute-error contract.\n\"\"\"\nfrom __future__ import annotations\n\nimport argparse\nimport csv\nimport gc\nimport hashlib\nimport itertools\nimport json\nimport math\nimport os\nimport platform\nimport re\nimport resource\nimport shutil\nimport subprocess\nimport sys\nimport tarfile\nimport tempfile\nimport time\nimport zipfile\nfrom dataclasses import dataclass, field\nfrom pathlib import Path\nfrom typing import Dict, Iterable, List, Optional, Sequence, Tuple\n\nimport numpy as np\nimport pandas as pd\n\n# -----------------------------------------------------------------------------\n# KAGGLE PASTE SETTINGS — edit only these two lines if needed.\n# First clean paper run: True. If Kaggle interrupts after scientific points have\n# been checkpointed, change to False before pasting/running the cell again.\n# -----------------------------------------------------------------------------\nKAGGLE_FRESH_FROM_ZERO = False\nKAGGLE_REQUIRE_BOTH_CODECS = True\n\n# -----------------------------------------------------------------------------\n# Configuration\n# -----------------------------------------------------------------------------\nROOT = Path(os.environ.get(\"HOC_WORK\", \"/kaggle/working/tbd_dimension_layout_final\"))\nOUT = ROOT / \"evidence\"\n# Software cache is deliberately outside the scientific run directory. A --fresh\n# run deletes every measurement/checkpoint but reuses binaries whose source commit\n# and executable are verified, avoiding minutes of unnecessary rebuild time.\nTOOLS = Path(os.environ.get(\"HOC_TOOLS\", \"/kaggle/working/high_order_compression_tools_cache\"))\nCACHE = ROOT / \"cache\"\nfor p in (ROOT, OUT, TOOLS, CACHE):\n p.mkdir(parents=True, exist_ok=True)\n\nINPUT_ROOT = Path(os.environ.get(\"HOC_INPUT_ROOT\", \"/kaggle/input\"))\nSEED = int(os.environ.get(\"HOC_SEED\", \"20260916\"))\nRNG = np.random.default_rng(SEED)\n\n# Headline quality grid. Each value is converted to an absolute tolerance from\n# the development-only normalized data range. Final comparisons always use the\n# measured reconstruction PSNR/RMSE, not the nominal value.\nNOMINAL_PSNR_GRID = np.asarray([60, 70, 80, 90, 100, 110, 120], dtype=float)\nSEARCH_NOMINAL_GRID = np.asarray([70, 100], dtype=float)\nVALID_NOMINAL_GRID = np.asarray([65, 80, 95, 110], dtype=float)\n\nPROXY_GLOBAL_KEEP = int(os.environ.get(\"HOC_PROXY_GLOBAL_KEEP\", \"24\"))\nPROXY_PER_PARTITION_KEEP = int(os.environ.get(\"HOC_PROXY_PER_PARTITION_KEEP\", \"48\"))\nCOARSE_KEEP = int(os.environ.get(\"HOC_COARSE_KEEP\", \"24\"))\nFINAL_KEEP = int(os.environ.get(\"HOC_FINAL_KEEP\", \"12\"))\nPROBE_MAX_VALUES = int(os.environ.get(\"HOC_PROBE_MAX_VALUES\", \"900000\"))\nCOMMON_WRAPPER_BYTES = int(os.environ.get(\"HOC_WRAPPER_BYTES\", \"64\"))\nSTRICT_DATASETS = os.environ.get(\"HOC_STRICT_DATASETS\", \"1\") != \"0\"\nSTRICT_CODECS = os.environ.get(\"HOC_STRICT_CODECS\", \"1\") != \"0\"\nMRI_CROP_ZYX = tuple(int(x) for x in os.environ.get(\"HOC_MRI_CROP_ZYX\", \"64,128,128\").split(\",\"))\nMRI_GROUP_SIZE = int(os.environ.get(\"HOC_MRI_GROUP_SIZE\", \"1\"))\nMAX_ODC_SUBJECTS = int(os.environ.get(\"HOC_MAX_ODC_SUBJECTS\", \"30\"))\nHYPERLEAF_MAX_PLOTS = int(os.environ.get(\"HOC_HYPERLEAF_MAX_PLOTS\", \"24\"))\nHYPERLEAF_LEAVES_PER_PLOT = int(os.environ.get(\"HOC_HYPERLEAF_LEAVES_PER_PLOT\", \"6\"))\nHYPERLEAF_PLOTS_PER_TEST_UNIT = int(os.environ.get(\"HOC_HYPERLEAF_PLOTS_PER_TEST_UNIT\", \"1\"))\nHYPERLEAF_CROP_YX = tuple(int(x) for x in os.environ.get(\"HOC_HYPERLEAF_CROP_YX\", \"32,128\").split(\",\"))\nIXI_MAX_SUBJECTS = int(os.environ.get(\"HOC_IXI_MAX_SUBJECTS\", \"24\"))\nIXI_GROUP_SIZE = int(os.environ.get(\"HOC_IXI_GROUP_SIZE\", \"1\"))\nMAX_DIAGNOSTIC_CHANNELS = int(os.environ.get(\"HOC_MAX_DIAGNOSTIC_CHANNELS\", \"16\"))\nBUILD_JOBS = max(1, int(os.environ.get(\"HOC_BUILD_JOBS\", \"2\")))\nGEOMETRY_MAX_VALUES = int(os.environ.get(\"HOC_GEOMETRY_MAX_VALUES\", \"500000\"))\nRESUME_ENABLED = os.environ.get(\"HOC_RESUME\", \"1\") != \"0\"\nif \"HOC_FRESH\" in os.environ:\n FRESH_RUN = os.environ.get(\"HOC_FRESH\", \"0\") == \"1\"\nelse:\n FRESH_RUN = bool(KAGGLE_FRESH_FROM_ZERO)\nFORCE_REBUILD = os.environ.get(\"HOC_FORCE_REBUILD\", \"0\") == \"1\"\nCHECKPOINT_SCHEMA = \"hoc-tbd-five-domain-uniform-selection-v2-stratified-intervals\"\nFINAL_CODECS = (\"SZ3\", \"HPEZ\")\nSELECTION_REVISION = \"uniform-stratified-measured-psnr-validation-v2\"\nVALID_WINNER_GRID = int(os.environ.get(\"HOC_VALID_WINNER_GRID\", \"41\"))\nAMS_TEST_WINDOW = int(os.environ.get(\"HOC_AMS_TEST_WINDOW\", \"256\"))\nAMS_MAX_TEST_WINDOWS = int(os.environ.get(\"HOC_AMS_MAX_TEST_WINDOWS\", \"6\"))\nBOUND_ACCEPT_RATIO = float(os.environ.get(\"HOC_BOUND_ACCEPT_RATIO\", \"1.0001\"))\n# Uniform descriptive measured-PSNR analysis. These are reporting constants, not layout-selection inputs.\nPSNR_INTERVAL_START_DB = float(os.environ.get(\"HOC_PSNR_INTERVAL_START_DB\", \"65\"))\nPSNR_INTERVAL_WIDTH_DB = float(os.environ.get(\"HOC_PSNR_INTERVAL_WIDTH_DB\", \"10\"))\nPSNR_INTERVAL_COUNT = int(os.environ.get(\"HOC_PSNR_INTERVAL_COUNT\", \"6\"))\nPSNR_INTERVAL_GRID_POINTS = int(os.environ.get(\"HOC_PSNR_INTERVAL_GRID_POINTS\", \"101\"))\nPSNR_INTERVAL_MIN_COVERAGE_DB = float(os.environ.get(\"HOC_PSNR_INTERVAL_MIN_COVERAGE_DB\", \"5\"))\n# AMS held-out evaluation uses the separate Kaggle competition GEFS test archive.\nAMS_CONFIRMATORY_WINDOW = int(os.environ.get(\"HOC_AMS_TEST_WINDOW\", \"256\"))\nAMS_CONFIRMATORY_MAX_WINDOWS = int(os.environ.get(\"HOC_AMS_TEST_WINDOWS\", \"6\"))\nAMS_CONFIRMATORY_SOURCE = \"gefs_test\"\n\nDATASET_URLS = {\n \"AMS/GEFS\": \"https://www.kaggle.com/competitions/ams-2014-solar-energy-prediction-contest/data\",\n \"Longitudinal MR-MS\": \"https://www.kaggle.com/datasets/farahmo/longitudinal-mri-data\",\n \"INSTANT-ODC multimodal MRI\": \"https://www.kaggle.com/competitions/instant-odc-ai-hackathon\",\n \"HyperLeaf2024 hyperspectral\": \"https://www.kaggle.com/competitions/HyperLeaf2024/data\",\n \"Preprocessed IXI MRI\": \"https://www.kaggle.com/datasets/hamedamin/preprocessed-oasis-and-epilepsy-and-ixi\",\n}\n\nHPEZ_REPO = \"https://github.com/szcompressor/QoZ.git\"\nHPEZ_COMMIT = os.environ.get(\"HOC_HPEZ_COMMIT\", \"a506cb0b1599f810c6f56e75edfdc418f9fd226f\")\n\nLOG_PATH = OUT / \"run_log.txt\"\nEVENTS_PATH = OUT / \"progress_events.csv\"\nSTATUS_PATH = OUT / \"run_status.json\"\nRESUME_STATE_PATH = OUT / \"resume_state.json\"\nFINGERPRINT_PATH = OUT / \"resume_protocol_fingerprint.json\"\nCHECKPOINT_TABLES = (\n \"dataset_units\", \"normalization\", \"layout_search\", \"layout_rankings\",\n \"selected_layouts\", \"raw_rd_points\", \"per_channel_metrics\",\n \"geometry_metrics\", \"codec_eligibility\",\n)\n\n# The extension executes only the two backends that showed useful representation sensitivity in the completed study.\nREQUIRED_CODECS = (\"SZ3\", \"HPEZ\")\nOPTIONAL_CODECS = set()\nSTRICT_ALL_CODEC_GATE = (os.environ.get(\"HOC_STRICT_ALL_CODEC_GATE\", \"1\") != \"0\") and bool(KAGGLE_REQUIRE_BOTH_CODECS)\nRUN_K_ABLATION = False\nPROXY_LAYOUT_CACHE: Dict[Tuple[str,int], List[Tuple[float, \"Layout\"]]] = {}\n\n\ndef log(msg: str) -> None:\n print(msg, flush=True)\n with LOG_PATH.open(\"a\", encoding=\"utf-8\") as f:\n f.write(msg + \"\\n\")\n\n\nRUN_T0 = time.perf_counter()\n\n\n\n\ndef _fmt_elapsed(seconds: float) -> str:\n seconds = max(0, int(seconds))\n h, rem = divmod(seconds, 3600)\n m, sec = divmod(rem, 60)\n return f\"{h:02d}:{m:02d}:{sec:02d}\"\n\ndef _memory_text() -> str:\n try:\n vals = {}\n for line in Path('/proc/meminfo').read_text().splitlines():\n if ':' in line:\n k, v = line.split(':', 1)\n vals[k] = int(v.strip().split()[0]) * 1024\n total = vals.get('MemTotal', 0) / 2**30\n avail = vals.get('MemAvailable', 0) / 2**30\n used = max(0.0, total-avail)\n return f\"RAM used={used:.1f} GiB, available={avail:.1f}/{total:.1f} GiB\"\n except Exception:\n return \"RAM unavailable\"\n\ndef progress_message(phase: str, action: str, detail: str = \"\") -> None:\n elapsed = _fmt_elapsed(time.perf_counter()-RUN_T0)\n msg = f\"[{phase}] {action} | elapsed={elapsed} | {_memory_text()}\"\n if detail:\n msg += f\"\\n {detail}\"\n log(msg)\n\ndef _tail_text(path: Path, n: int = 4) -> str:\n try:\n lines = path.read_text(errors='replace').splitlines()\n return \" | \".join(x.strip() for x in lines[-n:] if x.strip())[-1200:]\n except Exception:\n return \"\"\n\ndef run_command_progress(cmd: Sequence[str], label: str, cwd: Optional[Path] = None,\n timeout: Optional[int] = None, env: Optional[dict] = None,\n heartbeat: int = 15, check: bool = False) -> subprocess.CompletedProcess:\n \"\"\"Run a long external command with visible heartbeats instead of silent blocking.\"\"\"\n safe = re.sub(r'[^A-Za-z0-9_.-]+', '_', label)[:80]\n out = OUT / f\"command_{safe}.txt\"\n progress_message(\"SETUP\", f\"START {label}\", \"Command: \" + \" \".join(map(str, cmd)))\n t0 = time.perf_counter()\n with out.open('w', encoding='utf-8') as fh:\n proc = subprocess.Popen([str(x) for x in cmd], cwd=None if cwd is None else str(cwd),\n stdout=fh, stderr=subprocess.STDOUT, text=True, env=env)\n while proc.poll() is None:\n if timeout is not None and time.perf_counter()-t0 > timeout:\n proc.kill(); proc.wait()\n raise TimeoutError(f\"{label} exceeded {timeout} seconds; see {out.name}\")\n time.sleep(max(2, heartbeat))\n detail = _tail_text(out, 3) or \"No new command output yet; process is still running.\"\n progress_message(\"SETUP\", f\"RUNNING {label} ({_fmt_elapsed(time.perf_counter()-t0)})\", detail)\n rc = int(proc.returncode)\n text = out.read_text(errors='replace') if out.exists() else ''\n result = subprocess.CompletedProcess([str(x) for x in cmd], rc, text, None)\n if rc == 0:\n progress_message(\"SETUP\", f\"OK {label}\", f\"Finished in {_fmt_elapsed(time.perf_counter()-t0)}. Log: {out.name}\")\n else:\n progress_message(\"SETUP\", f\"FAIL {label}\", (_tail_text(out, 8) or f\"return code {rc}\"))\n if check:\n raise RuntimeError(f\"{label} failed (return code {rc}); see {out}\\n{text[-4000:]}\")\n return result\n\n\ndef emit_event(stage: str, status: str, dataset: str = \"\", compressor: str = \"\",\n unit: str = \"\", detail: str = \"\") -> None:\n \"\"\"Append a machine-readable progress event and refresh run_status.json.\"\"\"\n row = {\n \"timestamp_utc\": time.strftime(\"%Y-%m-%dT%H:%M:%SZ\", time.gmtime()),\n \"stage\": stage, \"status\": status, \"dataset\": dataset,\n \"compressor\": compressor, \"unit\": unit, \"detail\": detail[:1000],\n }\n exists = EVENTS_PATH.exists()\n with EVENTS_PATH.open(\"a\", newline=\"\", encoding=\"utf-8\") as f:\n w = csv.DictWriter(f, fieldnames=list(row))\n if not exists:\n w.writeheader()\n w.writerow(row)\n try:\n ev = pd.read_csv(EVENTS_PATH)\n counts = ev[\"status\"].value_counts().to_dict() if len(ev) else {}\n except Exception:\n counts = {}\n STATUS_PATH.write_text(json.dumps({\"last_event\": row, \"status_counts\": counts}, indent=2), encoding=\"utf-8\")\n\n\ndef _protocol_fingerprint_payload() -> dict:\n \"\"\"Fields that must match before scientific measurements can be reused.\"\"\"\n return {\n \"schema\": CHECKPOINT_SCHEMA,\n \"seed\": SEED,\n \"nominal_psnr_grid\": NOMINAL_PSNR_GRID.tolist(),\n \"search_nominal_grid\": SEARCH_NOMINAL_GRID.tolist(),\n \"valid_nominal_grid\": VALID_NOMINAL_GRID.tolist(),\n \"proxy_global_keep\": PROXY_GLOBAL_KEEP,\n \"proxy_per_partition_keep\": PROXY_PER_PARTITION_KEEP,\n \"coarse_keep\": COARSE_KEEP,\n \"final_keep\": FINAL_KEEP,\n \"probe_max_values\": PROBE_MAX_VALUES,\n \"common_wrapper_bytes\": COMMON_WRAPPER_BYTES,\n \"mri_crop_zyx\": list(MRI_CROP_ZYX),\n \"mri_group_size\": MRI_GROUP_SIZE,\n \"max_odc_subjects\": MAX_ODC_SUBJECTS,\n \"hyperleaf_max_plots\": HYPERLEAF_MAX_PLOTS,\n \"hyperleaf_leaves_per_plot\": HYPERLEAF_LEAVES_PER_PLOT,\n \"hyperleaf_plots_per_test_unit\": HYPERLEAF_PLOTS_PER_TEST_UNIT,\n \"hyperleaf_crop_yx\": list(HYPERLEAF_CROP_YX),\n \"ams_confirmatory_source\": AMS_CONFIRMATORY_SOURCE,\n \"ams_confirmatory_window\": AMS_CONFIRMATORY_WINDOW,\n \"ams_confirmatory_max_windows\": AMS_CONFIRMATORY_MAX_WINDOWS,\n \"selection_revision\": SELECTION_REVISION,\n \"final_codecs\": list(FINAL_CODECS),\n \"hpez_commit\": HPEZ_COMMIT,\n \"ixi_max_subjects\": IXI_MAX_SUBJECTS,\n \"ixi_group_size\": IXI_GROUP_SIZE,\n \"ams_test_window\": AMS_TEST_WINDOW,\n \"ams_max_test_windows\": AMS_MAX_TEST_WINDOWS,\n \"run_k_ablation\": RUN_K_ABLATION,\n \"psnr_interval_start_db\": PSNR_INTERVAL_START_DB,\n \"psnr_interval_width_db\": PSNR_INTERVAL_WIDTH_DB,\n \"psnr_interval_count\": PSNR_INTERVAL_COUNT,\n \"psnr_interval_grid_points\": PSNR_INTERVAL_GRID_POINTS,\n \"psnr_interval_min_coverage_db\": PSNR_INTERVAL_MIN_COVERAGE_DB,\n }\n\n\ndef protocol_fingerprint() -> str:\n blob = json.dumps(_protocol_fingerprint_payload(), sort_keys=True, separators=(\",\", \":\")).encode()\n return hashlib.sha256(blob).hexdigest()\n\n\ndef parse_layout_machine_label(text: str) -> Layout:\n m = re.fullmatch(r\"perm=([0-9,]+);cuts=([0-9,]*)\", str(text).strip())\n if not m:\n raise ValueError(f\"invalid layout label: {text!r}\")\n perm = tuple(int(x) for x in m.group(1).split(\",\") if x != \"\")\n cuts = tuple(int(x) for x in m.group(2).split(\",\") if x != \"\")\n return Layout(perm, cuts)\n\n\ndef _records_from_csv(path: Path) -> List[dict]:\n try:\n if not path.exists() or path.stat().st_size == 0:\n return []\n df = pd.read_csv(path)\n return df.where(pd.notna(df), None).to_dict(\"records\")\n except Exception as e:\n progress_message(\"RESUME AUDIT\", f\"Checkpoint ignored: {path.name}\", str(e))\n return []\n\n\ndef input_inventory_signature() -> str:\n \"\"\"Cheap identity for attached Kaggle inputs using relative paths and sizes only.\"\"\"\n h = hashlib.sha256()\n if not INPUT_ROOT.exists():\n h.update(b\"\")\n return h.hexdigest()\n count = 0\n for q in sorted((x for x in INPUT_ROOT.rglob(\"*\") if x.is_file()), key=lambda x: str(x)):\n try:\n rel = str(q.relative_to(INPUT_ROOT)).replace(os.sep, \"/\")\n st = q.stat()\n h.update(rel.encode(\"utf-8\", \"surrogatepass\")); h.update(b\"\\0\")\n h.update(str(int(st.st_size)).encode()); h.update(b\"\\n\")\n count += 1\n except Exception:\n continue\n h.update(f\"count={count}\".encode())\n return h.hexdigest()\n\n\ndef load_resume_state() -> Tuple[Dict[str, List[dict]], dict, bool]:\n \"\"\"Load checkpoints only when the exact scientific protocol matches.\"\"\"\n empty = {k: [] for k in CHECKPOINT_TABLES}\n inv_sig = input_inventory_signature()\n state = {\"fingerprint\": protocol_fingerprint(), \"input_inventory_signature\": inv_sig, \"completed_datasets\": {}}\n if not RESUME_ENABLED:\n progress_message(\"RESUME AUDIT\", \"Resume disabled\", \"HOC_RESUME=0; all scientific measurements will be recomputed.\")\n return empty, state, False\n current = protocol_fingerprint()\n previous = None\n try:\n if FINGERPRINT_PATH.exists():\n previous = json.loads(FINGERPRINT_PATH.read_text()).get(\"fingerprint\")\n except Exception:\n previous = None\n if previous != current:\n why = \"no compatible fingerprint exists\" if previous is None else \"stored fingerprint differs from current protocol\"\n stale_files=[OUT/f\"{name}.csv\" for name in CHECKPOINT_TABLES if (OUT/f\"{name}.csv\").exists()]\n if stale_files:\n stamp=time.strftime(\"%Y%m%d_%H%M%S\", time.localtime())\n stale_dir=OUT/f\"stale_checkpoint_{stamp}\"; stale_dir.mkdir(parents=True,exist_ok=True)\n for q in stale_files:\n try: shutil.move(str(q), str(stale_dir/q.name))\n except Exception: pass\n progress_message(\"RESUME AUDIT\", \"Incompatible scientific tables quarantined\", f\"Moved {len(stale_files)} table(s) to {stale_dir.name}; they cannot be mistaken for current results.\")\n progress_message(\"RESUME AUDIT\", \"Scientific checkpoints will NOT be reused\", f\"{why}. Existing tools/builds will still be verified and reused when valid.\")\n return empty, state, False\n loaded_state = None\n try:\n if RESUME_STATE_PATH.exists():\n candidate = json.loads(RESUME_STATE_PATH.read_text())\n if candidate.get(\"fingerprint\") == current:\n loaded_state = candidate\n except Exception as e:\n progress_message(\"RESUME AUDIT\", \"Resume-state JSON could not be read\", str(e))\n if loaded_state is None or loaded_state.get(\"input_inventory_signature\") != inv_sig:\n reason = \"resume_state.json is missing/incompatible\" if loaded_state is None else \"attached Kaggle input inventory changed\"\n stale_files=[OUT/f\"{name}.csv\" for name in CHECKPOINT_TABLES if (OUT/f\"{name}.csv\").exists()]\n if stale_files:\n stamp=time.strftime(\"%Y%m%d_%H%M%S\", time.localtime())\n stale_dir=OUT/f\"stale_input_checkpoint_{stamp}\"; stale_dir.mkdir(parents=True,exist_ok=True)\n for q in stale_files:\n try: shutil.move(str(q), str(stale_dir/q.name))\n except Exception: pass\n progress_message(\"RESUME AUDIT\", \"Scientific tables quarantined because inputs changed\", f\"Moved {len(stale_files)} table(s) to {stale_dir.name}.\")\n progress_message(\"RESUME AUDIT\", \"Scientific checkpoints will NOT be reused\", f\"{reason}. This prevents results from a different dataset attachment being mixed into the current run.\")\n return empty, state, False\n state = loaded_state\n for name in CHECKPOINT_TABLES:\n empty[name] = _records_from_csv(OUT / f\"{name}.csv\")\n counts = \", \".join(f\"{k}={len(v)}\" for k, v in empty.items() if len(v)) or \"no table rows\"\n progress_message(\"RESUME AUDIT\", \"Compatible scientific checkpoint found\", counts)\n return empty, state, True\n\n\ndef _dedupe_records(rows: List[dict], keys: Sequence[str]) -> List[dict]:\n out = {}\n for r in rows:\n key = tuple(str(r.get(k, \"\")) for k in keys)\n out[key] = r\n return list(out.values())\n\n\ndef checkpoint_tables(tables: Dict[str, List[dict]]) -> None:\n \"\"\"Persist enough state to resume after a Kaggle restart/failure.\"\"\"\n keymap = {\n \"dataset_units\": (\"dataset\",\"unit\",\"role\"),\n \"normalization\": (\"dataset\",\"channel\"),\n \"layout_search\": (\"dataset\",\"compressor\",\"stage\",\"layout\"),\n \"layout_rankings\": (\"dataset\",\"compressor\",\"rank\"),\n \"selected_layouts\": (\"dataset\",\"compressor\"),\n \"raw_rd_points\": (\"dataset\",\"unit\",\"compressor\",\"layout_kind\",\"nominal_psnr\"),\n \"per_channel_metrics\": (\"dataset\",\"unit\",\"compressor\",\"layout_kind\",\"nominal_psnr\",\"channel\"),\n \"geometry_metrics\": (\"dataset\",\"unit\",\"compressor\",\"layout_kind\",\"presented_axis\"),\n \"codec_eligibility\": (\"dataset\",\"compressor\"),\n }\n for name in CHECKPOINT_TABLES:\n rows = tables.get(name, [])\n if rows:\n rows = _dedupe_records(rows, keymap[name])\n tables[name] = rows\n pd.DataFrame(rows).to_csv(OUT / f\"{name}.csv\", index=False)\n payload = {\"fingerprint\": protocol_fingerprint(), \"payload\": _protocol_fingerprint_payload()}\n FINGERPRINT_PATH.write_text(json.dumps(payload, indent=2), encoding=\"utf-8\")\n\n\ndef save_resume_state(state: dict) -> None:\n state[\"fingerprint\"] = protocol_fingerprint()\n state.setdefault(\"input_inventory_signature\", input_inventory_signature())\n RESUME_STATE_PATH.write_text(json.dumps(state, indent=2), encoding=\"utf-8\")\n\n\ndef successful_rd_key_exists(tables: Dict[str, List[dict]], dataset: str, unit: str,\n compressor: str, kind: str, db: float,\n layout_label: str, version: str) -> bool:\n for r in tables.get(\"raw_rd_points\", []):\n try:\n if (str(r.get(\"dataset\")) == dataset and str(r.get(\"unit\")) == unit and\n str(r.get(\"compressor\")) == compressor and str(r.get(\"layout_kind\")) == kind and\n abs(float(r.get(\"nominal_psnr\")) - float(db)) < 1e-9 and\n str(r.get(\"layout\")) == layout_label and str(r.get(\"compressor_version\")) == str(version) and\n r.get(\"bits_per_value\") is not None and np.isfinite(float(r.get(\"bits_per_value\"))) and\n np.isfinite(float(r.get(\"bound_ratio\", np.inf))) and float(r.get(\"bound_ratio\", np.inf)) <= BOUND_ACCEPT_RATIO):\n return True\n except Exception:\n continue\n return False\n\n\ndef resumed_layout(bundle: DatasetBundle, adapter: Adapter,\n tables: Dict[str, List[dict]]) -> Optional[Tuple[Layout,List[Layout]]]:\n version = str(adapter.version())\n selected_row = None\n for r in tables.get(\"selected_layouts\", []):\n if str(r.get(\"dataset\")) == bundle.dataset_id and str(r.get(\"compressor\")) == adapter.name:\n if (str(r.get(\"compressor_version\", version)) == version and\n str(r.get(\"selection_revision\", \"\")) == SELECTION_REVISION):\n selected_row = r\n if selected_row is None:\n return None\n try:\n selected = parse_layout_machine_label(str(selected_row[\"selected_layout\"]))\n ranked_rows = [r for r in tables.get(\"layout_rankings\", [])\n if str(r.get(\"dataset\")) == bundle.dataset_id and str(r.get(\"compressor\")) == adapter.name\n and str(r.get(\"compressor_version\", version)) == version\n and str(r.get(\"selection_revision\", \"\")) == SELECTION_REVISION]\n ranked_rows.sort(key=lambda r: int(float(r.get(\"rank\", 999999))))\n ranked = [parse_layout_machine_label(str(r[\"layout\"])) for r in ranked_rows]\n if not ranked:\n ranked = [selected]\n progress_message(\"RESUME\", f\"Reusing frozen layout for {bundle.dataset_id}/{adapter.name}\",\n f\"selected={selected.label(bundle.axis_names)}; validated ranking rows={len(ranked)}\")\n return selected, ranked\n except Exception as e:\n progress_message(\"RESUME\", f\"Stored layout rejected for {bundle.dataset_id}/{adapter.name}\", str(e))\n return None\n\n\ndef _adapter_signature(adapters: Sequence[Adapter]) -> dict:\n return {a.name: str(a.version()) for a in adapters}\n\ndef dataset_marked_complete(state: dict, dataset_id: str, adapters: Sequence[Adapter]) -> bool:\n entry=state.get(\"completed_datasets\", {}).get(dataset_id, None)\n # Legacy boolean completion markers are intentionally treated as incomplete; per-point checkpoints are still reused.\n if not isinstance(entry, dict): return False\n return entry.get(\"compressors\", {}) == _adapter_signature(adapters)\n\ndef mark_dataset_complete(state: dict, dataset_id: str, adapters: Sequence[Adapter]) -> None:\n state.setdefault(\"completed_datasets\", {})[dataset_id] = {\"compressors\": _adapter_signature(adapters)}\n save_resume_state(state)\n\n\n\ndef prune_nonfinal_codec_rows(tables: Dict[str, List[dict]]) -> None:\n \"\"\"Remove stale rows from codecs excluded from this fresh four-codec experiment.\"\"\"\n allowed=set(FINAL_CODECS)\n for name, rows in list(tables.items()):\n if not isinstance(rows, list) or not rows:\n continue\n kept=[]\n for row in rows:\n if isinstance(row, dict) and \"compressor\" in row and str(row.get(\"compressor\")) not in allowed:\n continue\n kept.append(row)\n tables[name]=kept\n\ndef _plot_progress(tables: Dict[str, List[dict]], current_dataset: str = \"\") -> None:\n \"\"\"Write diagnostic PNGs from the measurements accumulated so far.\"\"\"\n try:\n import matplotlib\n matplotlib.use(\"Agg\")\n import matplotlib.pyplot as plt\n except Exception as e:\n log(f\"PLOT WARN: matplotlib unavailable: {e}\")\n return\n\n raw = pd.DataFrame(tables.get(\"raw_rd_points\", []))\n if len(raw) and \"compressor\" in raw.columns:\n tmp = raw.copy()\n bpv = tmp[\"bits_per_value\"] if \"bits_per_value\" in tmp else pd.Series(np.nan, index=tmp.index)\n tmp[\"result\"] = np.where(bpv.notna(), \"successful\", \"failed\")\n ct = tmp.groupby([\"compressor\", \"result\"]).size().unstack(fill_value=0)\n for col in (\"successful\", \"failed\"):\n if col not in ct:\n ct[col] = 0\n fig, ax = plt.subplots(figsize=(7.0, 4.2))\n x = np.arange(len(ct.index)); ok = ct[\"successful\"].to_numpy(); bad = ct[\"failed\"].to_numpy()\n ax.bar(x, ok, label=\"successful\")\n ax.bar(x, bad, bottom=ok, label=\"failed\")\n ax.set_xticks(x, ct.index, rotation=25, ha=\"right\")\n ax.set_ylabel(\"Measured RD points\")\n ax.set_title(\"Benchmark progress by compressor\")\n ax.legend()\n fig.tight_layout(); fig.savefig(OUT/\"progress_completion.png\", dpi=140); plt.close(fig)\n\n needed = {\"dataset\",\"unit\",\"compressor\",\"layout_kind\",\"psnr_db\",\"bits_per_value\"}\n if len(raw) and needed.issubset(raw.columns):\n q = raw[raw[\"bits_per_value\"].notna() & raw[\"psnr_db\"].notna()].copy()\n if len(q):\n ds = current_dataset if current_dataset in set(q.dataset) else str(q.iloc[-1].dataset)\n q = q[q.dataset == ds]\n unit = str(q.iloc[-1].unit); q = q[q.unit == unit]\n fig, ax = plt.subplots(figsize=(7.2, 4.6))\n for (codec, kind), g in q.groupby([\"compressor\",\"layout_kind\"]):\n g = g.sort_values(\"psnr_db\")\n ax.plot(g.psnr_db, g.bits_per_value, marker=\"o\", label=f\"{codec} {kind}\")\n ax.set_xlabel(\"Measured PSNR (dB)\"); ax.set_ylabel(\"Bits/value\")\n ax.set_title(f\"Current rate-distortion curves: {ds} / {unit}\")\n ax.grid(True, alpha=0.25); ax.legend(fontsize=7, ncol=2)\n fig.tight_layout(); fig.savefig(OUT/\"progress_rd_current.png\", dpi=140); plt.close(fig)\n\n if len(raw) and {\"dataset\",\"compressor\",\"bound_ratio\"}.issubset(raw.columns):\n q = raw[raw[\"bound_ratio\"].notna()].copy()\n if len(q):\n mx = q.groupby([\"dataset\",\"compressor\"], as_index=False).bound_ratio.max()\n labels = [f\"{a}/{b}\" for a,b in zip(mx.dataset, mx.compressor)]\n fig, ax = plt.subplots(figsize=(max(7.0, .55*len(mx)), 4.2))\n ax.bar(np.arange(len(mx)), mx.bound_ratio.to_numpy())\n ax.axhline(1.0, linestyle=\"--\", linewidth=1.2, label=\"requested bound\")\n ax.set_xticks(np.arange(len(mx)), labels, rotation=45, ha=\"right\")\n ax.set_ylabel(r\"max $\\|e\\|_\\infty / \\varepsilon$\")\n ax.set_title(\"Accumulated pointwise error-bound check\")\n ax.legend()\n fig.tight_layout(); fig.savefig(OUT/\"progress_error_bound.png\", dpi=140); plt.close(fig)\n\n\ndef _plot_final_summaries() -> None:\n \"\"\"Create compact final diagnostic plots after postprocessing.\"\"\"\n try:\n import matplotlib\n matplotlib.use(\"Agg\")\n import matplotlib.pyplot as plt\n except Exception:\n return\n p = OUT/\"rate_saving_summary.csv\"\n if p.exists():\n df = pd.read_csv(p)\n if len(df):\n labels = [f\"{d}/{c}\" for d,c in zip(df.dataset, df.compressor)]\n vals = df.mean_saving_pct.to_numpy(float)\n lo = vals - df.bootstrap95_lo.to_numpy(float); hi = df.bootstrap95_hi.to_numpy(float) - vals\n fig, ax = plt.subplots(figsize=(max(7.0, .55*len(df)), 4.5))\n ax.bar(np.arange(len(df)), vals, yerr=np.vstack([np.maximum(lo,0),np.maximum(hi,0)]), capsize=3)\n ax.axhline(0.0, linewidth=1.0)\n ax.set_xticks(np.arange(len(df)), labels, rotation=45, ha=\"right\")\n ax.set_ylabel(\"Integrated rate saving (%)\")\n ax.set_title(\"Validation-selected representation vs reference layout\")\n fig.tight_layout(); fig.savefig(OUT/\"final_rate_savings.png\", dpi=160); plt.close(fig)\n p = OUT/\"runtime_summary.csv\"\n if p.exists():\n df = pd.read_csv(p)\n if len(df) and \"layout_kind\" in df:\n df = df[df.layout_kind == \"selected\"]\n if len(df):\n targets = df.execution_target.astype(str) if \"execution_target\" in df.columns else pd.Series([\"unknown\"]*len(df))\n labels = [f\"{d}/{c} [{t}]\" for d,c,t in zip(df.dataset, df.compressor, targets)]\n fig, ax = plt.subplots(figsize=(max(7.0,.62*len(df)),4.7))\n ax.bar(np.arange(len(df)), df.median_encode_MBps.to_numpy(float))\n ax.set_xticks(np.arange(len(df)), labels, rotation=45, ha=\"right\")\n ax.set_ylabel(\"Median encode throughput (MB/s)\")\n ax.set_title(\"Selected-layout encode throughput\")\n fig.tight_layout(); fig.savefig(OUT/\"final_encode_throughput.png\", dpi=160); plt.close(fig)\n\n\ndef sh(cmd: Sequence[str], cwd: Optional[Path] = None, timeout: Optional[int] = None,\n env: Optional[dict] = None, check: bool = False) -> subprocess.CompletedProcess:\n p = subprocess.run([str(x) for x in cmd], cwd=None if cwd is None else str(cwd),\n stdout=subprocess.PIPE, stderr=subprocess.STDOUT,\n text=True, timeout=timeout, env=env)\n if check and p.returncode != 0:\n raise RuntimeError(\"command failed: \" + \" \".join(map(str, cmd)) + \"\\n\" + p.stdout[-4000:])\n return p\n\n\ndef package_version(name: str) -> str:\n try:\n import importlib.metadata\n return importlib.metadata.version(name)\n except Exception:\n return \"unavailable\"\n\n\ndef _unique_existing_dirs(paths: Sequence[Path]) -> List[Path]:\n out=[]; seen=set()\n for q in paths:\n try:\n q=Path(q).resolve()\n except Exception:\n q=Path(q)\n key=str(q)\n if key in seen or not q.is_dir():\n continue\n seen.add(key); out.append(q)\n return out\n\n\n\n\n\n# -----------------------------------------------------------------------------\n# Layouts\n# -----------------------------------------------------------------------------\n@dataclass(frozen=True)\nclass Layout:\n perm: Tuple[int, ...]\n cuts: Tuple[int, ...]\n\n @property\n def d(self) -> int:\n return len(self.perm)\n\n @property\n def m(self) -> int:\n return len(self.cuts) + 1\n\n @property\n def groups(self) -> Tuple[Tuple[int, int], ...]:\n b = (0,) + self.cuts + (self.d,)\n return tuple((b[i], b[i + 1]) for i in range(len(b) - 1))\n\n def label(self, axis_names: Sequence[str]) -> str:\n parts = []\n for lo, hi in self.groups:\n parts.append(\"\".join(axis_names[self.perm[j]] for j in range(lo, hi)))\n return \"|\".join(parts)\n\n def machine_label(self) -> str:\n return \"perm=\" + \",\".join(map(str, self.perm)) + \";cuts=\" + \",\".join(map(str, self.cuts))\n\n\ndef enumerate_layouts(d: int, m: int) -> List[Layout]:\n if not (1 <= m <= d):\n raise ValueError((d, m))\n if m == d:\n cuts = tuple(range(1, d))\n return [Layout(tuple(p), cuts) for p in itertools.permutations(range(d))]\n cutsets = list(itertools.combinations(range(1, d), m - 1))\n return [Layout(tuple(p), tuple(c)) for p in itertools.permutations(range(d)) for c in cutsets]\n\n\ndef reference_layout(d: int, m: int) -> Layout:\n # Native semantic order. The first d-m+1 axes are folded together; remaining\n # axes stay explicit. This is deterministic and contains no data-dependent choice.\n if m == d:\n cuts = tuple(range(1, d))\n else:\n first = d - m + 1\n cuts = tuple(range(first, d))\n return Layout(tuple(range(d)), cuts)\n\n\ndef folded_shape(shape: Sequence[int], layout: Layout) -> Tuple[int, ...]:\n perm_shape = [int(shape[i]) for i in layout.perm]\n out = []\n for lo, hi in layout.groups:\n out.append(int(np.prod(perm_shape[lo:hi], dtype=np.int64)))\n return tuple(out)\n\n\ndef fold_tensor(x: np.ndarray, layout: Layout) -> np.ndarray:\n y = np.transpose(x, layout.perm)\n return np.ascontiguousarray(y.reshape(folded_shape(x.shape, layout)), dtype=np.float32)\n\n\ndef unfold_tensor(y: np.ndarray, original_shape: Sequence[int], layout: Layout) -> np.ndarray:\n pshape = tuple(int(original_shape[i]) for i in layout.perm)\n z = np.asarray(y, dtype=np.float32).reshape(pshape)\n inv = tuple(np.argsort(np.asarray(layout.perm)))\n # Return a semantic-order view. Measurement code supports non-contiguous arrays,\n # so there is no reason to duplicate the entire reconstruction here.\n return np.transpose(z, inv)\n\n\ndef center_crop_probe(x: np.ndarray, max_values: int = PROBE_MAX_VALUES) -> np.ndarray:\n if x.size <= max_values:\n return np.ascontiguousarray(x, dtype=np.float32)\n original = list(x.shape)\n target = list(x.shape)\n while int(np.prod(target, dtype=np.int64)) > max_values:\n ax = int(np.argmax(target))\n if target[ax] <= 2:\n break\n target[ax] = max(2, int(math.ceil(target[ax] * 0.75)))\n slices = []\n for n, t in zip(original, target):\n start = max(0, (n - t) // 2)\n slices.append(slice(start, start + t))\n return np.ascontiguousarray(x[tuple(slices)], dtype=np.float32)\n\n\ndef proxy_tensor(x: np.ndarray, max_per_axis: int = 12) -> np.ndarray:\n sl = []\n for n in x.shape:\n if n <= max_per_axis:\n sl.append(slice(None))\n else:\n step = max(1, n // max_per_axis)\n sl.append(slice(0, n, step))\n return np.ascontiguousarray(x[tuple(sl)], dtype=np.float32)\n\n\ndef _locality_score_folded(y: np.ndarray, scale: float) -> float:\n \"\"\"Exact locality score used by the proxy stage for an already folded view.\n\n `scale` is invariant under permutations/foldings because every layout contains\n exactly the same proxy samples. Keeping it outside the inner loop avoids\n recomputing a full-array standard deviation thousands of times.\n \"\"\"\n vals = []\n for ax, n in enumerate(y.shape):\n if n > 1:\n vals.append(float(np.mean(np.abs(np.diff(y, axis=ax)), dtype=np.float64)) / scale)\n return float(np.mean(vals)) if vals else float(\"inf\")\n\n\ndef locality_score(x: np.ndarray, layout: Layout) -> float:\n # Retained for tests/small one-off calls. The full search uses\n # `proxy_layout_scores_fast`, which is mathematically identical but avoids\n # reconstructing the same proxy and permutation for every cut pattern.\n px = proxy_tensor(x)\n y = fold_tensor(px, layout)\n scale = float(np.std(px, dtype=np.float64)) + 1e-12\n return _locality_score_folded(y, scale)\n\n\ndef proxy_layout_scores_fast(x: np.ndarray, layouts: Sequence[Layout],\n axis_names: Sequence[str], dataset_id: str,\n compressor: str) -> List[Tuple[float, Layout]]:\n \"\"\"Evaluate the *same* proxy score for every layout with ~10x less copying.\n\n For d=6,m=4 there are 720 permutations and 10 cut patterns per\n permutation. The old implementation rebuilt/transposed a contiguous proxy\n 7200 times. Here each permutation is materialized once and its cut patterns\n are reshape views. No candidate or score definition is changed.\n \"\"\"\n if not layouts:\n return []\n px = proxy_tensor(x)\n scale = float(np.std(px, dtype=np.float64)) + 1e-12\n by_perm: Dict[Tuple[int, ...], List[Layout]] = {}\n for lay in layouts:\n by_perm.setdefault(tuple(lay.perm), []).append(lay)\n perms = sorted(by_perm)\n total = len(layouts)\n done = 0\n t0 = time.perf_counter()\n rows: List[Tuple[float, Layout]] = []\n report_every = max(1, int(os.environ.get(\"HOC_PROXY_REPORT_EVERY_PERMS\", \"25\")))\n top_show = max(1, int(os.environ.get(\"HOC_PROXY_TOP_SHOW\", \"5\")))\n log(f\"[PROXY SEARCH] {dataset_id}/{compressor}: evaluating {total} layouts \"\n f\"from {len(perms)} unique permutations on proxy shape={tuple(px.shape)}\")\n for pi, perm in enumerate(perms, 1):\n # One copy per permutation; all foldings below are views of this buffer.\n yp = np.ascontiguousarray(np.transpose(px, perm), dtype=np.float32)\n pshape = tuple(int(v) for v in yp.shape)\n for lay in by_perm[perm]:\n out_shape = []\n for lo, hi in lay.groups:\n out_shape.append(int(np.prod(pshape[lo:hi], dtype=np.int64)))\n yf = yp.reshape(tuple(out_shape))\n sc = _locality_score_folded(yf, scale)\n rows.append((sc, lay))\n done += 1\n if pi == 1 or pi % report_every == 0 or pi == len(perms):\n elapsed = max(time.perf_counter() - t0, 1e-9)\n rate = done / elapsed\n pct = 100.0 * done / max(total, 1)\n eta = (total - done) / rate if rate > 0 else float(\"nan\")\n best = sorted(rows, key=lambda z: (z[0], z[1].machine_label()))[:top_show]\n best_txt = \"; \".join(\n f\"{lay.label(axis_names)}={sc:.5f}\" for sc, lay in best\n )\n log(f\"[PROXY PROGRESS] {dataset_id}/{compressor}: {done}/{total} \"\n f\"({pct:.1f}%) | {rate:.1f} layouts/s | elapsed={elapsed:.1f}s | \"\n f\"ETA~{eta:.1f}s\")\n log(f\" current best {top_show}: {best_txt}\")\n rows.sort(key=lambda z: (z[0], z[1].machine_label()))\n elapsed = time.perf_counter() - t0\n log(f\"[PROXY DONE] {dataset_id}/{compressor}: {total} layouts in {elapsed:.1f}s; \"\n f\"best={rows[0][1].label(axis_names)} score={rows[0][0]:.6f}\")\n return rows\n\n# -----------------------------------------------------------------------------\n# Dataset structures and normalization\n# -----------------------------------------------------------------------------\n@dataclass\nclass TestUnit:\n unit_id: str\n x: np.ndarray\n roi_mask: Optional[np.ndarray] = None\n\n@dataclass\nclass DatasetBundle:\n dataset_id: str\n dataset_name: str\n axis_names: Tuple[str, ...]\n channel_axis: int\n channel_names: Tuple[str, ...]\n search: np.ndarray\n validation: np.ndarray\n tests: List[TestUnit]\n normalization_mean: np.ndarray\n normalization_std: np.ndarray\n normalization_ignore_zero: bool\n source_detail: str\n native_unit_label: str = \"native\"\n\n @property\n def d(self) -> int:\n return self.search.ndim\n\n\ndef channel_stats(x: np.ndarray, channel_axis: int, ignore_zero: bool) -> Tuple[np.ndarray, np.ndarray]:\n means, stds = [], []\n for c in range(x.shape[channel_axis]):\n sl = [slice(None)] * x.ndim\n sl[channel_axis] = c\n a = np.asarray(x[tuple(sl)], dtype=np.float64)\n m = np.isfinite(a)\n if ignore_zero:\n m &= np.abs(a) > 1e-12\n if not np.any(m):\n mu, sd = 0.0, 1.0\n else:\n mu = float(np.mean(a[m]))\n sd = max(float(np.std(a[m])), 1e-12)\n means.append(mu); stds.append(sd)\n return np.asarray(means), np.asarray(stds)\n\n\ndef normalize_channels(x: np.ndarray, channel_axis: int, mu: np.ndarray, sd: np.ndarray) -> np.ndarray:\n # One full float32 copy, then normalize channel-by-channel in place. This avoids\n # the several full-size temporaries created by a vectorized subtract/divide chain.\n out = np.array(x, dtype=np.float32, order=\"C\", copy=True)\n for c in range(out.shape[channel_axis]):\n sl = [slice(None)] * out.ndim\n sl[channel_axis] = c\n v = out[tuple(sl)]\n v -= np.float32(mu[c])\n v /= np.float32(sd[c])\n return out\n\n\ndef development_range(bundle: DatasetBundle) -> float:\n lo = min(float(np.nanmin(bundle.search)), float(np.nanmin(bundle.validation)))\n hi = max(float(np.nanmax(bundle.search)), float(np.nanmax(bundle.validation)))\n return max(hi - lo, 1e-12)\n\n\ndef eps_from_nominal(db: float, dr: float) -> float:\n return float(dr * 10.0 ** (-float(db) / 20.0))\n\n# -----------------------------------------------------------------------------\n# AMS / GEFS loader\n# -----------------------------------------------------------------------------\nAMS_VARS = (\"dswrf_sfc\", \"dlwrf_sfc\", \"tmp_2m\", \"spfh_2m\", \"pres_msl\", \"tcdc_eatm\")\n\ndef _find_ams_archive(split: str = \"train\") -> Path:\n split = str(split).strip().lower()\n if split not in {\"train\", \"test\"}:\n raise ValueError(f\"AMS split must be train or test, got {split!r}\")\n names = [f\"gefs_{split}.zip\", f\"gefs_{split}.tar.gz\"]\n cands = []\n for name in names:\n cands.extend(INPUT_ROOT.rglob(name))\n if not cands:\n raise FileNotFoundError(\n f\"AMS/GEFS {split} archive missing. Attach the Kaggle AMS 2013-2014 Solar Energy Prediction Contest data; expected one of {names}.\"\n )\n return sorted(cands, key=lambda p: len(str(p)))[0]\n\n\ndef _extract_ams_vars(archive: Path, split: str) -> Dict[str, Path]:\n # Train/test must have separate extraction directories because member basenames overlap.\n outdir = CACHE / f\"ams_netcdf_{split}\"\n outdir.mkdir(parents=True, exist_ok=True)\n result: Dict[str, Path] = {}\n for v in AMS_VARS:\n existing = list(outdir.glob(v + \"*.nc\"))\n if existing:\n result[v] = existing[0]\n if len(result) == len(AMS_VARS):\n return result\n if archive.suffix == \".zip\":\n with zipfile.ZipFile(archive) as zf:\n names = zf.namelist()\n for v in AMS_VARS:\n if v in result:\n continue\n matches = [n for n in names if Path(n).name.startswith(v + \"_\") and n.endswith(\".nc\")]\n if not matches:\n raise FileNotFoundError(f\"{v} not found in {archive}\")\n member = sorted(matches)[0]\n dst = outdir / Path(member).name\n with zf.open(member) as src, dst.open(\"wb\") as fh:\n shutil.copyfileobj(src, fh)\n result[v] = dst\n else:\n with tarfile.open(archive, \"r:gz\") as tf:\n members = [m for m in tf.getmembers() if m.isfile()]\n for v in AMS_VARS:\n if v in result:\n continue\n matches = [m for m in members if Path(m.name).name.startswith(v + \"_\") and m.name.endswith(\".nc\")]\n if not matches:\n raise FileNotFoundError(f\"{v} not found in {archive}\")\n m = sorted(matches, key=lambda z: z.name)[0]\n dst = outdir / Path(m.name).name\n src = tf.extractfile(m)\n if src is None:\n raise RuntimeError(m.name)\n with src, dst.open(\"wb\") as fh:\n shutil.copyfileobj(src, fh)\n result[v] = dst\n return result\n\n\ndef _ams_date_count(path: Path) -> int:\n from netCDF4 import Dataset\n with Dataset(path, \"r\") as ds:\n candidates = []\n for name, var in ds.variables.items():\n if getattr(var, \"ndim\", 0) == 5 and np.issubdtype(np.dtype(var.dtype), np.number):\n candidates.append((int(np.prod(var.shape)), name, tuple(int(n) for n in var.shape)))\n if not candidates:\n raise ValueError(f\"no numeric 5D field in {path}\")\n _, _, shape = max(candidates, key=lambda z: z[0])\n fixed = {11, 5, 9, 16}\n date_axes = [n for n in shape if n not in fixed]\n if len(date_axes) != 1:\n raise ValueError(f\"cannot infer AMS date axis from {shape} in {path}\")\n return int(date_axes[0])\n\ndef _read_ams_slice(path: Path, start: int, count: int) -> np.ndarray:\n from netCDF4 import Dataset\n with Dataset(path, \"r\") as ds:\n candidates = []\n for name, var in ds.variables.items():\n if getattr(var, \"ndim\", 0) == 5 and np.issubdtype(np.dtype(var.dtype), np.number):\n candidates.append((int(np.prod(var.shape)), name, var))\n if not candidates: raise ValueError(f\"no numeric 5D field in {path}\")\n _, _, var = max(candidates, key=lambda z: z[0])\n shape = tuple(int(n) for n in var.shape)\n def axis_of(size: int) -> int:\n ids = [i for i, n in enumerate(shape) if n == size]\n if len(ids) != 1: raise ValueError((shape, size))\n return ids[0]\n ens, lead, lat, lon = axis_of(11), axis_of(5), axis_of(9), axis_of(16)\n date = [i for i in range(5) if i not in {ens, lead, lat, lon}][0]\n sl = [slice(None)] * 5; sl[date] = slice(start, start + count)\n raw = var[tuple(sl)]\n if np.ma.isMaskedArray(raw): raw = np.ma.filled(raw, np.nan)\n raw = np.asarray(raw, dtype=np.float32)\n raw = np.transpose(raw, axes=(date, ens, lead, lat, lon))\n return raw\n\n\ndef load_ams() -> DatasetBundle:\n train_archive = _find_ams_archive(\"train\")\n test_archive = _find_ams_archive(\"test\")\n train_sources = _extract_ams_vars(train_archive, \"train\")\n test_sources = _extract_ams_vars(test_archive, \"test\")\n\n # Development and validation use the competition training archive; held-out evaluation uses the separate competition-test archive.\n norm_range = (0, 128)\n search_range = (128, 512)\n valid_range = (640, 512)\n\n stats = []\n for v in AMS_VARS:\n a = _read_ams_slice(train_sources[v], *norm_range)\n finite = np.isfinite(a)\n stats.append((float(np.mean(a[finite])), max(float(np.std(a[finite])), 1e-12)))\n mu = np.asarray([x[0] for x in stats])\n sd = np.asarray([x[1] for x in stats])\n\n def tensor(sources: Dict[str, Path], rng: Tuple[int, int]) -> np.ndarray:\n chans = []\n for i, v in enumerate(AMS_VARS):\n a = _read_ams_slice(sources[v], *rng)\n if a.shape[0] != rng[1]:\n raise RuntimeError(f\"AMS slice {rng} for {v} returned {a.shape[0]} dates\")\n a = np.where(np.isfinite(a), a, mu[i]).astype(np.float32)\n chans.append((a - np.float32(mu[i])) / np.float32(sd[i]))\n return np.ascontiguousarray(np.stack(chans, axis=0), dtype=np.float32)\n\n search = tensor(train_sources, search_range)\n valid = tensor(train_sources, valid_range)\n\n counts = [_ams_date_count(test_sources[v]) for v in AMS_VARS]\n total_test_dates = min(counts)\n if len(set(counts)) != 1:\n log(f\"WARNING AMS test variables have unequal date counts {dict(zip(AMS_VARS, counts))}; using common minimum={total_test_dates}\")\n win = AMS_CONFIRMATORY_WINDOW\n nwin = min(AMS_CONFIRMATORY_MAX_WINDOWS, total_test_dates // win)\n if nwin < 1:\n raise RuntimeError(f\"AMS confirmatory test archive has only {total_test_dates} dates; need at least {win}\")\n # Spread non-overlapping confirmatory windows across the full independent test archive.\n # No quality/layout decision may depend on these windows. When the archive has spare\n # dates, the spacing leaves gaps between blocks, reducing short-range temporal coupling.\n if nwin == 1:\n starts = [0]\n else:\n starts = np.linspace(0, total_test_dates - win, nwin).round().astype(int).tolist()\n # Guard against rounding overlap for unusually short archives.\n for i in range(1, len(starts)):\n starts[i] = max(starts[i], starts[i-1] + win)\n if starts[-1] + win > total_test_dates:\n shift = starts[-1] + win - total_test_dates\n starts = [max(0, int(s - shift)) for s in starts]\n test_ranges = [(int(s), win) for s in starts]\n tests = [TestUnit(f\"competition_test_{s}_{s+n}\", tensor(test_sources, (s, n))) for s, n in test_ranges]\n progress_message(\n \"AMS SPLIT\",\n \"Independent competition-test confirmation configured\",\n f\"development=train[{search_range[0]}:{search_range[0]+search_range[1]}], validation=train[{valid_range[0]}:{valid_range[0]+valid_range[1]}], \"\n f\"confirmatory source={test_archive.name}, available_dates={total_test_dates}, windows={test_ranges}; all held-out endpoints use the separate competition-test archive\"\n )\n source_detail = f\"development={train_archive}; confirmatory={test_archive}\"\n return DatasetBundle(\"ams_gefs\", \"AMS/GEFS\", (\"V\",\"D\",\"E\",\"F\",\"Y\",\"X\"), 0,\n AMS_VARS, search, valid, tests, mu, sd, False, source_detail, \"physical\")\n\n# -----------------------------------------------------------------------------\n# MRI utilities and loaders\n# -----------------------------------------------------------------------------\ndef _load_nii_zyx(path: Path) -> np.ndarray:\n import nibabel as nib\n a = np.asarray(nib.load(str(path)).get_fdata(dtype=np.float32), dtype=np.float32)\n if a.ndim != 3: raise ValueError(f\"expected 3D NIfTI, got {a.shape}: {path}\")\n # NIfTI array is conventionally X,Y,Z in numpy; present spatial axes as Z,Y,X.\n return np.ascontiguousarray(np.transpose(a, (2,1,0)), dtype=np.float32)\n\n\ndef crop_or_pad_zyx(a: np.ndarray, target: Tuple[int,int,int], fill: float = 0.0) -> np.ndarray:\n out = np.full(target, fill, dtype=a.dtype)\n src_sl, dst_sl = [], []\n for n, t in zip(a.shape, target):\n take = min(n, t)\n s0 = max(0, (n - take)//2); d0 = max(0, (t - take)//2)\n src_sl.append(slice(s0, s0+take)); dst_sl.append(slice(d0, d0+take))\n out[tuple(dst_sl)] = a[tuple(src_sl)]\n return out\n\n\ndef deterministic_split(ids: List[str], max_items: Optional[int] = None) -> Tuple[List[str],List[str],List[str]]:\n ids = sorted(ids)\n rng = np.random.default_rng(SEED)\n perm = list(np.asarray(ids, dtype=object)[rng.permutation(len(ids))])\n if max_items is not None: perm = perm[:max_items]\n n = len(perm)\n if n < 8: raise RuntimeError(f\"need >=8 subjects/patients, found {n}\")\n n_search = max(3, int(round(0.40*n)))\n n_val = max(2, int(round(0.20*n)))\n if n_search + n_val > n - 2:\n n_val = max(2, n - n_search - 2)\n return perm[:n_search], perm[n_search:n_search+n_val], perm[n_search+n_val:]\n\n\ndef group_ids(ids: List[str], group_size: int) -> List[List[str]]:\n return [ids[i:i+group_size] for i in range(0, len(ids), group_size) if len(ids[i:i+group_size]) == group_size]\n\n\ndef _mrms_nifti_files(root: Path) -> List[Path]:\n \"\"\"Return NIfTI files below *root* without assuming Kaggle's mount basename.\"\"\"\n out: List[Path] = []\n try:\n for p in root.rglob(\"*\"):\n if not p.is_file():\n continue\n n = p.name.lower()\n if n.endswith(\".nii\") or n.endswith(\".nii.gz\"):\n out.append(p)\n except Exception:\n pass\n return out\n\n\ndef _mrms_candidate_roots() -> List[Path]:\n \"\"\"Locate plausible Kaggle mounts for the Longitudinal MR_MS dataset.\"\"\"\n if not INPUT_ROOT.exists():\n return []\n roots = [p for p in INPUT_ROOT.iterdir() if p.is_dir()]\n keys = (\"longitudinal\", \"mr-ms\", \"mr_ms\", \"longitudinal-mri\", \"farah\", \"multiple-sclerosis\")\n preferred = [p for p in roots if any(k in p.name.lower() for k in keys)]\n # The canonical Kaggle slug currently mounts as longitudinal-mri-data, but\n # Kaggle may normalize source names differently. Prefer plausible roots,\n # then fall back to every attached source so discovery is robust.\n ordered: List[Path] = []\n seen = set()\n for p in preferred + roots:\n rp = str(p.resolve())\n if rp not in seen:\n ordered.append(p); seen.add(rp)\n return ordered\n\n\ndef _extract_mrms_archives_if_needed(roots: Sequence[Path]) -> List[Path]:\n \"\"\"Extract a dataset archive only when no mounted NIfTI files are visible.\"\"\"\n archive_candidates: List[Path] = []\n for root in roots:\n try:\n for p in root.rglob(\"*.zip\"):\n low = p.name.lower()\n whole = str(p).lower()\n if any(k in whole for k in (\"longitudinal\", \"mr_ms\", \"mr-ms\", \"mri\", \"patient\")):\n archive_candidates.append(p)\n except Exception:\n pass\n if not archive_candidates:\n return []\n dst = CACHE / \"mrms_extracted\"\n marker = dst / \".complete\"\n if marker.exists():\n return [dst]\n dst.mkdir(parents=True, exist_ok=True)\n for z in archive_candidates:\n try:\n if zipfile.is_zipfile(z):\n print(f\"[MR-MS DISCOVERY] Mounted NIfTI files were not visible; extracting archive: {z}\", flush=True)\n with zipfile.ZipFile(z) as zh:\n zh.extractall(dst)\n except Exception as e:\n print(f\"[MR-MS DISCOVERY] Archive extraction failed for {z}: {e}\", flush=True)\n if _mrms_nifti_files(dst):\n marker.write_text(\"ok\\n\")\n return [dst]\n return []\n\n\ndef _find_mrms_patients() -> Dict[str, Dict[str, Path]]:\n \"\"\"Discover complete two-visit T1/T2/FLAIR patients robustly.\n\n Accepts both .nii and .nii.gz, case differences, nested Kaggle mount\n directories, optional ``_reg`` suffixes, and a raw/registered copy. When\n both raw and registered images exist, registered non-raw images are\n preferred. This does not alter the scientific split or tensor definition.\n \"\"\"\n roots = _mrms_candidate_roots()\n top = [p.name for p in roots[:20]]\n print(f\"[MR-MS DISCOVERY] Searching Kaggle inputs. Candidate roots={top}\", flush=True)\n\n nifti: List[Path] = []\n used_roots: List[Path] = []\n # Search plausible roots first and stop once the intended dataset is found;\n # this avoids traversing the entire BraTS tree unnecessarily.\n for root in roots:\n fs = _mrms_nifti_files(root)\n if fs:\n # Identify the source by the actual longitudinal filename grammar,\n # not by the mount-directory title.\n longitudinal = [p for p in fs if re.search(r\"study[_-]?[12]\", p.name, re.I)]\n multimodal = [p for p in longitudinal if re.search(r\"(?:T1|T2|FLAIR)\", p.name, re.I)]\n if multimodal:\n nifti.extend(fs); used_roots.append(root)\n # One Kaggle source is enough; do not accidentally combine MRI datasets.\n break\n\n if not nifti:\n extra = _extract_mrms_archives_if_needed(roots)\n for root in extra:\n fs = _mrms_nifti_files(root)\n if fs:\n nifti.extend(fs); used_roots.append(root)\n\n if used_roots:\n print(f\"[MR-MS DISCOVERY] Detected source root: {used_roots[0]}\", flush=True)\n print(f\"[MR-MS DISCOVERY] NIfTI files visible in detected source: {len(nifti)}\", flush=True)\n if nifti:\n sample = \", \".join(str(p.relative_to(used_roots[0])) if used_roots and used_roots[0] in p.parents else p.name for p in nifti[:8])\n print(f\"[MR-MS DISCOVERY] Example files: {sample}\", flush=True)\n\n # Flexible grammar for the published dataset convention:\n # patientX_study1_T1_reg.nii.gz, etc. ``reg`` is optional so a Kaggle\n # version containing only raw files can still be diagnosed/used.\n img_re = re.compile(\n r\"^(?P.+?)[_-]study[_-]?(?P[12])[_-](?PT1|T2|FLAIR)(?P.*?)(?:\\\\.nii(?:\\\\.gz)?)$\",\n re.I,\n )\n # The raw regex above is intentionally reconstructed below after stripping\n # extensions; this avoids ambiguity around '.nii.gz'.\n core_re = re.compile(r\"^(?P.+?)[_-]study[_-]?(?P[12])[_-](?PT1|T2|FLAIR)(?P.*)$\", re.I)\n\n candidates: Dict[str, Dict[str, Tuple[int, Path]]] = {}\n masks: Dict[str, Dict[str, Path]] = {}\n for p in nifti:\n name = p.name\n low = name.lower()\n core = re.sub(r\"\\.nii(?:\\.gz)?$\", \"\", name, flags=re.I)\n m = core_re.match(core)\n if m:\n pid = m.group(\"pid\")\n visit = int(m.group(\"visit\"))\n mod = m.group(\"mod\").upper()\n tail = m.group(\"tail\").lower()\n key = f\"s{visit}_{mod}\"\n # Registered, non-raw images are the intended published inputs.\n score = 0\n if \"reg\" in tail: score += 10\n if \"raw\" not in [q.lower() for q in p.parts]: score += 5\n old = candidates.setdefault(pid, {}).get(key)\n if old is None or score > old[0]:\n candidates[pid][key] = (score, p)\n continue\n # Optional ROI/brain masks.\n mm = re.match(r\"^(?P.+?)[_-](?Pbrainmask|gt)$\", core, re.I)\n if mm:\n pid = mm.group(\"pid\"); kind = mm.group(\"kind\").lower()\n masks.setdefault(pid, {})[kind] = p\n\n patients: Dict[str, Dict[str, Path]] = {}\n required = [f\"s{v}_{m}\" for v in (1,2) for m in (\"T1\",\"T2\",\"FLAIR\")]\n for pid, d in candidates.items():\n if all(k in d for k in required):\n req: Dict[str, Path] = {k: d[k][1] for k in required}\n req.update(masks.get(pid, {}))\n patients[pid] = req\n\n print(f\"[MR-MS DISCOVERY] Complete patients with 2 visits x 3 modalities: {len(patients)}\", flush=True)\n if patients:\n first = sorted(patients)[0]\n print(f\"[MR-MS DISCOVERY] First complete patient={first}; files={len(patients[first])}\", flush=True)\n else:\n # Detailed diagnostic rather than the previous opaque FileNotFoundError.\n attached = []\n if INPUT_ROOT.exists():\n attached = sorted(p.name for p in INPUT_ROOT.iterdir())\n modality_counts = {}\n for pid, d in candidates.items():\n modality_counts[pid] = sorted(d)\n examples = list(modality_counts.items())[:5]\n raise FileNotFoundError(\n \"Longitudinal MR_MS dataset was not discoverable as complete two-visit \"\n \"T1/T2/FLAIR patients. The correct Kaggle source is \"\n \"farahmo/longitudinal-mri-data. \"\n f\"Attached /kaggle/input roots={attached[:30]}; \"\n f\"candidate NIfTI count={len(nifti)}; partial patient examples={examples}. \"\n \"If candidate NIfTI count is 0, the Kaggle dataset is not actually mounted \"\n \"in this notebook session even if its data-card page is open.\"\n )\n return patients\n\n\ndef load_mrms() -> DatasetBundle:\n pats = _find_mrms_patients()\n if not pats: raise FileNotFoundError(\"Longitudinal MR_MS dataset not found\")\n search_ids, val_ids, test_ids = deterministic_split(list(pats))\n mods = (\"T1\",\"T2\",\"FLAIR\")\n def patient(pid: str) -> Tuple[np.ndarray, Optional[np.ndarray]]:\n d = pats[pid]; visits=[]\n for visit in (1,2):\n vv=[]\n for mod in mods:\n vv.append(crop_or_pad_zyx(_load_nii_zyx(d[f\"s{visit}_{mod}\"]), MRI_CROP_ZYX))\n visits.append(np.stack(vv, axis=0))\n x = np.stack(visits, axis=0) # visit, modality, z,y,x\n mask = None\n if \"gt\" in d:\n mask = crop_or_pad_zyx(_load_nii_zyx(d[\"gt\"]), MRI_CROP_ZYX) > 0\n return x.astype(np.float32), mask\n def stack(ids: List[str]) -> np.ndarray:\n return np.stack([patient(i)[0] for i in ids], axis=0)\n raw_search = stack(search_ids); raw_val = stack(val_ids)\n mu, sd = channel_stats(raw_search, 2, True)\n search = normalize_channels(raw_search, 2, mu, sd); valid = normalize_channels(raw_val, 2, mu, sd)\n del raw_search, raw_val\n gc.collect()\n tests=[]\n for g in group_ids(test_ids, MRI_GROUP_SIZE):\n xs=[]; masks=[]\n for pid in g:\n x,m=patient(pid); xs.append(x); masks.append(m)\n x=np.stack(xs,axis=0); x=normalize_channels(x,2,mu,sd)\n mask = None if any(m is None for m in masks) else np.stack(masks,axis=0)\n tests.append(TestUnit(\"patients_\"+\"_\".join(g),x,mask))\n if not tests: raise RuntimeError(\"MR-MS split produced no full test group\")\n return DatasetBundle(\"mrms\",\"Longitudinal MR-MS\",(\"P\",\"S\",\"M\",\"Z\",\"Y\",\"X\"),2,mods,\n search,valid,tests,mu,sd,True,str(next(iter(pats.values()))[\"s1_T1\"]),\"MRI intensity\")\n\n\ndef _brats_subjects() -> Dict[str, Dict[str, Path]]:\n all_nii = list(INPUT_ROOT.rglob(\"*.nii\")) + list(INPUT_ROOT.rglob(\"*.nii.gz\"))\n bydir: Dict[Path, List[Path]] = {}\n for p in all_nii:\n bydir.setdefault(p.parent, []).append(p)\n out: Dict[str, Dict[str, Path]] = {}\n for d, files in bydir.items():\n low = {p.name.lower():p for p in files}\n def pick(pred):\n c=[p for p in files if pred(p.name.lower())]\n return sorted(c)[0] if c else None\n flair=pick(lambda n:\"flair\" in n and \"seg\" not in n)\n t2=pick(lambda n: re.search(r\"(^|_)t2($|_|\\.)\", n) is not None and \"flair\" not in n)\n t1ce=pick(lambda n: (\"t1ce\" in n or \"t1gd\" in n or \"t1c\" in n) and \"seg\" not in n)\n t1=pick(lambda n: re.search(r\"(^|_)t1($|_|\\.)\", n) is not None and not any(z in n for z in (\"t1ce\",\"t1gd\",\"t1c\",\"seg\")))\n if all(z is not None for z in (t1,t1ce,t2,flair)):\n seg=pick(lambda n:\"seg\" in n)\n pid=d.name\n out[pid]={\"T1\":t1,\"T1ce\":t1ce,\"T2\":t2,\"FLAIR\":flair}\n if seg is not None: out[pid][\"seg\"]=seg\n return out\n\n\ndef load_brats() -> DatasetBundle:\n subs=_brats_subjects()\n if not subs: raise FileNotFoundError(\"BraTS 2020/2021 multimodal NIfTI data not found\")\n search_ids,val_ids,test_ids=deterministic_split(list(subs),MAX_BRATS_SUBJECTS)\n mods=(\"T1\",\"T1ce\",\"T2\",\"FLAIR\")\n def subject(pid:str)->Tuple[np.ndarray,Optional[np.ndarray]]:\n d=subs[pid]; arr=[]\n for mod in mods: arr.append(crop_or_pad_zyx(_load_nii_zyx(d[mod]),MRI_CROP_ZYX))\n x=np.stack(arr,axis=0).astype(np.float32)\n mask=None\n if \"seg\" in d: mask=crop_or_pad_zyx(_load_nii_zyx(d[\"seg\"]),MRI_CROP_ZYX)>0\n return x,mask\n def stack(ids:List[str])->np.ndarray: return np.stack([subject(i)[0] for i in ids],axis=0)\n raw_search=stack(search_ids); raw_val=stack(val_ids)\n mu,sd=channel_stats(raw_search,1,True)\n search=normalize_channels(raw_search,1,mu,sd); valid=normalize_channels(raw_val,1,mu,sd)\n del raw_search, raw_val\n gc.collect()\n tests=[]\n for g in group_ids(test_ids,MRI_GROUP_SIZE):\n xs=[]; ms=[]\n for pid in g:\n x,m=subject(pid); xs.append(x); ms.append(m)\n x=normalize_channels(np.stack(xs,axis=0),1,mu,sd)\n mask=None if any(m is None for m in ms) else np.stack(ms,axis=0)\n tests.append(TestUnit(\"subjects_\"+\"_\".join(g),x,mask))\n if not tests: raise RuntimeError(\"BraTS split produced no full test group\")\n return DatasetBundle(\"brats\",\"BraTS multimodal MRI\",(\"P\",\"M\",\"Z\",\"Y\",\"X\"),1,mods,\n search,valid,tests,mu,sd,True,str(next(iter(subs.values()))[\"T1\"]),\"MRI intensity\")\n\n# -----------------------------------------------------------------------------\n# Cross-domain extension loaders\n# -----------------------------------------------------------------------------\n\ndef _strip_nii_ext(name: str) -> str:\n return re.sub(r\"\\.nii(?:\\.gz)?$\", \"\", str(name), flags=re.I)\n\n\ndef _instant_odc_subjects() -> Dict[str, Dict[str, Path]]:\n \"\"\"Discover one INSTANT-ODC/BraTS-style four-modality NIfTI source.\n\n The loader intentionally selects a single top-level Kaggle source rather than\n combining NIfTI files from unrelated attached datasets.\n \"\"\"\n if not INPUT_ROOT.exists():\n return {}\n roots=[p for p in INPUT_ROOT.iterdir() if p.is_dir()]\n def root_priority(r: Path) -> Tuple[int,str]:\n n=r.name.lower()\n score=0\n if \"instant\" in n: score+=100\n if \"odc\" in n: score+=80\n if \"hackathon\" in n: score+=30\n if \"brats\" in n: score+=10\n return (-score,n)\n roots=sorted(roots,key=root_priority)\n\n def classify(core: str) -> Tuple[Optional[str], Optional[str]]:\n # Return patient id and canonical modality. Longest/specific tokens first.\n pats=[\n (\"T1ce\", r\"(?:^|[_-])(t1ce|t1c|t1gd)(?:$|[_-])\"),\n (\"FLAIR\",r\"(?:^|[_-])(flair|t2f)(?:$|[_-])\"),\n (\"T2\", r\"(?:^|[_-])(t2|t2w)(?:$|[_-])\"),\n (\"T1\", r\"(?:^|[_-])(t1|t1n)(?:$|[_-])\"),\n (\"seg\", r\"(?:^|[_-])(seg|mask|label)(?:$|[_-])\"),\n ]\n low=core.lower()\n for mod,pat in pats:\n m=re.search(pat,low,re.I)\n if not m: continue\n # Remove exactly the modality token while preserving the rest of the id.\n a,b=m.span(1)\n pid=(core[:a]+core[b:]).strip(\"_- .\")\n pid=re.sub(r\"[_-]{2,}\",\"_\",pid)\n if pid:\n return pid,mod\n return None,None\n\n best: Dict[str, Dict[str, Path]]={}\n best_root=None\n for root in roots:\n files=[]\n try:\n files=[q for q in root.rglob(\"*\") if q.is_file() and (q.name.lower().endswith(\".nii\") or q.name.lower().endswith(\".nii.gz\"))]\n except Exception:\n continue\n if not files: continue\n subj: Dict[str, Dict[str, Path]]={}\n for q in files:\n core=_strip_nii_ext(q.name)\n pid,mod=classify(core)\n if pid is None or mod is None: continue\n subj.setdefault(pid,{})[mod]=q\n complete={pid:d for pid,d in subj.items() if all(m in d for m in (\"T1\",\"T1ce\",\"T2\",\"FLAIR\"))}\n print(f\"[INSTANT-ODC DISCOVERY] root={root.name}: NIfTI={len(files)}, complete four-modality subjects={len(complete)}\",flush=True)\n if len(complete)>len(best):\n best,best_root=complete,root\n # A strongly named source with enough subjects is preferred immediately.\n if len(complete)>=8 and any(k in root.name.lower() for k in (\"instant\",\"odc\")):\n best,best_root=complete,root\n break\n if best_root is not None:\n print(f\"[INSTANT-ODC DISCOVERY] selected root={best_root}; subjects={len(best)}\",flush=True)\n first=sorted(best)[0] if best else None\n if first:\n print(f\"[INSTANT-ODC DISCOVERY] example subject={first}; files={[best[first][m].name for m in ('T1','T1ce','T2','FLAIR')]}\",flush=True)\n return best\n\n\ndef load_instant_odc() -> DatasetBundle:\n subs=_instant_odc_subjects()\n if not subs:\n raise FileNotFoundError(\n \"INSTANT-ODC four-modality NIfTI data not found. Attach the Kaggle competition \"\n \"'INSTANT-ODC AI Hackathon' to this notebook.\"\n )\n search_ids,val_ids,test_ids=deterministic_split(list(subs),MAX_ODC_SUBJECTS)\n mods=(\"T1\",\"T1ce\",\"T2\",\"FLAIR\")\n def subject(pid:str)->Tuple[np.ndarray,Optional[np.ndarray]]:\n d=subs[pid]; arr=[]\n for mod in mods:\n arr.append(crop_or_pad_zyx(_load_nii_zyx(d[mod]),MRI_CROP_ZYX))\n x=np.stack(arr,axis=0).astype(np.float32)\n mask=None\n if \"seg\" in d:\n mask=crop_or_pad_zyx(_load_nii_zyx(d[\"seg\"]),MRI_CROP_ZYX)>0\n return x,mask\n def stack(ids:List[str])->np.ndarray:\n return np.stack([subject(i)[0] for i in ids],axis=0)\n raw_search=stack(search_ids); raw_val=stack(val_ids)\n mu,sd=channel_stats(raw_search,1,True)\n search=normalize_channels(raw_search,1,mu,sd)\n valid=normalize_channels(raw_val,1,mu,sd)\n del raw_search,raw_val; gc.collect()\n tests=[]\n for g in group_ids(test_ids,MRI_GROUP_SIZE):\n xs=[]; ms=[]\n for pid in g:\n x,m=subject(pid); xs.append(x); ms.append(m)\n x=normalize_channels(np.stack(xs,axis=0),1,mu,sd)\n mask=None if any(m is None for m in ms) else np.stack(ms,axis=0)\n tests.append(TestUnit(\"subjects_\"+\"_\".join(g),x,mask))\n if not tests:\n raise RuntimeError(\"INSTANT-ODC split produced no full held-out subject group\")\n src=str(next(iter(subs.values()))[\"T1\"])\n progress_message(\"INSTANT-ODC SPLIT\",\"Subject-disjoint extension configured\",\n f\"development subjects={len(search_ids)}; validation={len(val_ids)}; held-out subjects={len(test_ids)} in {len(tests)} groups; crop={MRI_CROP_ZYX}\")\n return DatasetBundle(\"instant_odc\",\"INSTANT-ODC multimodal MRI\",(\"P\",\"M\",\"Z\",\"Y\",\"X\"),1,mods,\n search,valid,tests,mu,sd,True,src,\"MRI intensity\")\n\n\ndef _center_crop_or_pad_yx(a: np.ndarray, target: Tuple[int,int]) -> np.ndarray:\n \"\"\"Center crop/pad the final two dimensions of BxYxX data.\"\"\"\n if a.ndim!=3:\n raise ValueError(f\"expected BxYxX hyperspectral cube, got {a.shape}\")\n b,y,x=a.shape; ty,tx=target\n out=np.zeros((b,ty,tx),dtype=a.dtype)\n sy=min(y,ty); sx=min(x,tx)\n ys=max(0,(y-sy)//2); xs=max(0,(x-sx)//2)\n yd=max(0,(ty-sy)//2); xd=max(0,(tx-sx)//2)\n out[:,yd:yd+sy,xd:xd+sx]=a[:,ys:ys+sy,xs:xs+sx]\n return out\n\n\ndef _hyperleaf_image_files(root: Path) -> List[Path]:\n \"\"\"Return HyperLeaf TIFF cubes without relying on extension case.\n\n The Kaggle competition stores train and test cubes together under an images\n directory. Some mounted copies preserve upper-case .TIF/.TIFF suffixes, for\n which pathlib.rglob(\"*.tiff\") is case-sensitive on Linux.\n \"\"\"\n files=[]\n for q in root.rglob(\"*\"):\n try:\n if q.is_file() and q.suffix.lower() in {\".tif\", \".tiff\"}:\n files.append(q)\n except OSError:\n continue\n return sorted(files, key=lambda q: str(q).lower())\n\n\ndef _find_case_insensitive(root: Path, basename: str) -> List[Path]:\n target=basename.lower()\n out=[]\n for q in root.rglob(\"*\"):\n try:\n if q.is_file() and q.name.lower()==target:\n out.append(q)\n except OSError:\n continue\n return sorted(out, key=lambda q: len(str(q)))\n\n\ndef _find_hyperleaf_root() -> Path:\n if not INPUT_ROOT.exists():\n raise FileNotFoundError(\"/kaggle/input missing\")\n roots=[p for p in INPUT_ROOT.iterdir() if p.is_dir()]\n preferred=[r for r in roots if \"hyperleaf\" in r.name.lower()]\n diagnostics=[]\n for root in preferred+[r for r in roots if r not in preferred]:\n trains=_find_case_insensitive(root,\"train.csv\")\n if not trains:\n continue\n tifs=_hyperleaf_image_files(root)\n diagnostics.append(f\"{root.name}: train.csv={len(trains)}, TIFF={len(tifs)}\")\n # Published competition has 1,590 train leaves and 820 test leaves. Use a\n # conservative threshold only to disambiguate it from unrelated train.csv files.\n if len(tifs)>=100:\n images_dirs=sorted({str(q.parent) for q in tifs[:200]})\n parent_hint=images_dirs[0] if images_dirs else \"(unknown)\"\n print(f\"[HYPERLEAF DISCOVERY] selected root={root}; train.csv={trains[0]}; TIFF files={len(tifs)}; example image dir={parent_hint}\",flush=True)\n return root\n detail=\" | \".join(diagnostics) if diagnostics else \"no root with train.csv was found\"\n raise FileNotFoundError(\n \"HyperLeaf2024 data not found/resolved under /kaggle/input. \"\n \"Expected the Kaggle competition mount containing train.csv and images/*.tif[f]. \"\n f\"Discovery: {detail}\")\n\n\ndef _read_hyperleaf_cube(path: Path) -> np.ndarray:\n import tifffile\n a=np.asarray(tifffile.imread(str(path)))\n a=np.squeeze(a)\n if a.ndim!=3:\n raise ValueError(f\"HyperLeaf TIFF must be 3-D, got {a.shape} at {path}\")\n # Published leaf cubes contain 204 spectral channels. Detect their axis rather than\n # assuming TIFF page/channel order.\n band_axes=[i for i,n in enumerate(a.shape) if int(n)==204]\n if not band_axes:\n # Tolerate transformed copies while insisting on a plausible hyperspectral axis.\n plausible=[i for i,n in enumerate(a.shape) if 64 <= int(n) <= 512]\n if not plausible:\n raise ValueError(f\"cannot identify spectral axis in HyperLeaf cube {a.shape}\")\n band_axis=min(plausible,key=lambda i:abs(int(a.shape[i])-204))\n else:\n band_axis=band_axes[0]\n a=np.moveaxis(a,band_axis,0) # B,Y,X\n if a.shape[1] > a.shape[2] and a.shape[1] > 256 and a.shape[2] <= 128:\n # Harmless orientation normalization for TIFF variants that swap the two spatial ranks.\n a=a.transpose(0,2,1)\n if np.issubdtype(a.dtype,np.integer):\n info=np.iinfo(a.dtype); scale=float(info.max) if info.max else 1.0\n a=a.astype(np.float32)/np.float32(scale)\n else:\n a=a.astype(np.float32,copy=False)\n a=np.where(np.isfinite(a),a,0.0).astype(np.float32,copy=False)\n return _center_crop_or_pad_yx(a,HYPERLEAF_CROP_YX)\n\n\ndef _hyperleaf_plot_groups(root: Path) -> Tuple[Dict[str,List[Path]],Path]:\n train_candidates=_find_case_insensitive(root,\"train.csv\")\n if not train_candidates:\n raise FileNotFoundError(\"HyperLeaf train.csv not found\")\n train_csv=train_candidates[0]\n df=pd.read_csv(train_csv)\n if \"ImageId\" not in df.columns:\n # Case-insensitive recovery.\n matches=[c for c in df.columns if str(c).strip().lower()==\"imageid\"]\n if not matches: raise ValueError(f\"HyperLeaf train.csv lacks ImageId; columns={list(df.columns)}\")\n df=df.rename(columns={matches[0]:\"ImageId\"})\n\n image_files=_hyperleaf_image_files(root)\n fmap={}\n numeric_fmap={}\n for q in image_files:\n # Kaggle ImageId is numeric in train.csv while the TIFF stems are zero-padded\n # (for example ImageId=0 maps to 00000.tiff). Keep both literal and canonical\n # integer indices so matching never depends on a particular padding width.\n for key in {q.name, q.stem, q.name.lower(), q.stem.lower()}:\n fmap[str(key).strip()]=q\n stem=q.stem.strip()\n try:\n # Accept stems such as 00000 and, defensively, 00000.0.\n sv=float(stem)\n if sv.is_integer():\n numeric_fmap[int(sv)]=q\n except Exception:\n pass\n\n def resolve(v) -> Optional[Path]:\n z=str(v).strip()\n # First try literal filename/stem forms.\n candidates=[z, z.lower(), Path(z).name, Path(z).name.lower(), Path(z).stem, Path(z).stem.lower()]\n for c in candidates:\n if c in fmap:\n return fmap[c]\n\n # Canonical numeric match: 0 -> 00000.tiff, 123 -> 00123.tiff,\n # 123.0 -> 00123.tiff, independent of the filename padding width.\n for raw in (z, Path(z).stem):\n try:\n fv=float(raw)\n if fv.is_integer():\n hit=numeric_fmap.get(int(fv))\n if hit is not None:\n return hit\n except Exception:\n pass\n\n # Last-resort suffix normalization.\n base=Path(z).stem.lower()\n return fmap.get(base)\n\n resolved=[resolve(v) for v in df.ImageId]\n df[\"_path\"]=resolved\n n_resolved=int(pd.Series(resolved).notna().sum())\n if n_resolved < 100:\n examples=[str(v) for v in df.ImageId.head(5).tolist()]\n file_examples=[q.name for q in image_files[:8]]\n raise RuntimeError(\n f\"HyperLeaf image resolution found only {n_resolved}/{len(df)} training TIFFs; \"\n f\"visible TIFFs={len(image_files)}; example ImageId={examples}; example files={file_examples}\")\n df=df[df[\"_path\"].notna()].copy()\n print(f\"[HYPERLEAF DISCOVERY] resolved {n_resolved}/{len(df)} train ImageIds against {len(image_files)} visible TIFF cubes\",flush=True)\n\n # Use an explicit plot column when a Kaggle revision exposes one. Otherwise recover\n # the published plot grouping from plot-level targets. GrainWeight, fertilizer, and\n # cultivar are constant at plot level in the acquisition design; Gsw/PhiPS2 are not\n # used for grouping.\n explicit=[c for c in df.columns if re.search(r\"(^|_)plot(id)?($|_)\",str(c),re.I)]\n if explicit:\n keys=df[explicit[0]].astype(str)\n grouping_mode=f\"explicit column {explicit[0]}\"\n else:\n req=[\"GrainWeight\",\"Fertilizer\",\"Heerup\",\"Kvium\",\"Rembrandt\",\"Sheriff\"]\n missing=[c for c in req if c not in df.columns]\n if missing:\n raise ValueError(f\"HyperLeaf plot grouping needs columns {req}; missing={missing}\")\n # Try several stable roundings and choose the grouping closest to the published 24 train plots.\n options=[]\n for dec in (8,6,5,4,3):\n vals=[]\n for _,r in df.iterrows():\n cult=max((\"Heerup\",\"Kvium\",\"Rembrandt\",\"Sheriff\"),key=lambda c:float(r[c]))\n vals.append((round(float(r[\"GrainWeight\"]),dec),round(float(r[\"Fertilizer\"]),dec),cult))\n n=len(set(vals)); options.append((abs(n-24),dec,vals,n))\n _,dec,vals,n=min(options,key=lambda z:(z[0],-z[1]))\n if not (16 <= n <= 40):\n raise RuntimeError(f\"could not reconstruct a plausible HyperLeaf plot partition; inferred groups={n}\")\n keys=pd.Series([f\"GW{v[0]}_F{v[1]}_{v[2]}\" for v in vals],index=df.index)\n grouping_mode=f\"plot-level target signature, GrainWeight/Fertilizer rounded to {dec} decimals\"\n groups: Dict[str,List[Path]]={}\n for key,q in zip(keys,df[\"_path\"]):\n groups.setdefault(str(key),[]).append(Path(q))\n groups={k:sorted(v,key=lambda q:q.name) for k,v in groups.items() if len(v)>=HYPERLEAF_LEAVES_PER_PLOT}\n sizes=sorted(len(v) for v in groups.values())\n print(f\"[HYPERLEAF DISCOVERY] grouping={grouping_mode}; usable plots={len(groups)}; leaves/plot min/median/max={sizes[0] if sizes else 0}/{int(np.median(sizes)) if sizes else 0}/{sizes[-1] if sizes else 0}\",flush=True)\n return groups,train_csv\n\n\ndef load_hyperleaf() -> DatasetBundle:\n root=_find_hyperleaf_root(); plots,train_csv=_hyperleaf_plot_groups(root)\n if len(plots)<8:\n raise RuntimeError(f\"HyperLeaf requires >=8 usable plot groups, found {len(plots)}\")\n search_ids,val_ids,test_ids=deterministic_split(list(plots),HYPERLEAF_MAX_PLOTS)\n L=HYPERLEAF_LEAVES_PER_PLOT\n def plot_tensor(pid:str)->np.ndarray:\n fs=plots[pid]\n # Evenly spaced deterministic selection avoids depending on an arbitrary prefix of a plot.\n ix=np.linspace(0,len(fs)-1,L).round().astype(int)\n cubes=[_read_hyperleaf_cube(fs[int(i)]) for i in ix]\n b=min(c.shape[0] for c in cubes)\n if b<16: raise RuntimeError(f\"implausibly few HyperLeaf bands for plot {pid}: {b}\")\n cubes=[c[:b] for c in cubes]\n return np.stack(cubes,axis=0).astype(np.float32) # leaf,band,y,x\n # Confirm one spectral size first; all selected plot tensors must agree.\n def stack_plot_ids(ids:List[str])->np.ndarray:\n arr=[plot_tensor(i) for i in ids]\n b=min(x.shape[1] for x in arr)\n return np.stack([x[:,:b] for x in arr],axis=0) # plot,leaf,band,y,x\n raw_search=stack_plot_ids(search_ids); raw_val=stack_plot_ids(val_ids)\n b=min(raw_search.shape[2],raw_val.shape[2])\n raw_search=raw_search[:,:,:b]; raw_val=raw_val[:,:,:b]\n mu,sd=channel_stats(raw_search,2,True)\n search=normalize_channels(raw_search,2,mu,sd); valid=normalize_channels(raw_val,2,mu,sd)\n del raw_search,raw_val; gc.collect()\n tests=[]\n for g in group_ids(test_ids,HYPERLEAF_PLOTS_PER_TEST_UNIT):\n x=stack_plot_ids(g)[:,:,:b]\n x=normalize_channels(x,2,mu,sd)\n tests.append(TestUnit(\"plots_\"+\"__\".join(g),x,None))\n if not tests:\n raise RuntimeError(\"HyperLeaf split produced no full held-out plot group\")\n channel_names=tuple(f\"band_{i:03d}\" for i in range(b))\n progress_message(\"HYPERLEAF SPLIT\",\"Plot-disjoint hyperspectral extension configured\",\n f\"development plots={len(search_ids)}; validation={len(val_ids)}; held-out plots={len(test_ids)} in {len(tests)} groups; leaves/plot={L}; bands={b}; crop_yx={HYPERLEAF_CROP_YX}\")\n return DatasetBundle(\"hyperleaf\",\"HyperLeaf2024 hyperspectral leaves\",(\"P\",\"L\",\"B\",\"Y\",\"X\"),2,channel_names,\n search,valid,tests,mu,sd,True,str(train_csv),\"reflectance\")\n\n\n\n# -----------------------------------------------------------------------------\n# Preprocessed IXI MRI loader\n# -----------------------------------------------------------------------------\n\ndef _ixi_candidate_roots() -> List[Path]:\n if not INPUT_ROOT.exists(): return []\n roots=[p for p in INPUT_ROOT.iterdir() if p.is_dir()]\n keys=(\"preprocessed\",\"ixi\",\"oasis\",\"epilepsy\",\"hamedamin\")\n preferred=[p for p in roots if any(k in p.name.lower() for k in keys)]\n return preferred+[p for p in roots if p not in preferred]\n\n\ndef _maybe_extract_ixi(root: Path) -> Path:\n nii=list(root.rglob(\"*.nii\"))+list(root.rglob(\"*.nii.gz\"))\n if len(nii)>=12: return root\n archives=list(root.rglob(\"*.zip\"))+list(root.rglob(\"*.tar\"))+list(root.rglob(\"*.tar.gz\"))+list(root.rglob(\"*.tgz\"))\n if not archives: return root\n dest=CACHE/\"ixi_extracted\"/root.name\n if dest.exists() and len(list(dest.rglob(\"*.nii\"))+list(dest.rglob(\"*.nii.gz\")))>=12:\n return dest\n dest.mkdir(parents=True,exist_ok=True)\n # Prefer archives whose names mention IXI; extract one archive at a time until NIfTI appears.\n archives=sorted(archives,key=lambda p:(0 if \"ixi\" in p.name.lower() else 1,p.name))\n for arc in archives[:4]:\n progress_message(\"IXI DISCOVERY\",f\"Extracting candidate archive {arc.name}\",f\"destination={dest}\")\n try:\n if arc.name.lower().endswith(\".zip\"):\n with zipfile.ZipFile(arc) as zf: zf.extractall(dest)\n else:\n with tarfile.open(arc) as tf: tf.extractall(dest)\n except Exception as e:\n log(f\"[IXI DISCOVERY] archive extraction skipped for {arc}: {type(e).__name__}: {e}\")\n if len(list(dest.rglob(\"*.nii\"))+list(dest.rglob(\"*.nii.gz\")))>=12: return dest\n return dest if dest.exists() else root\n\n\ndef _ixi_map_subjects() -> Tuple[Dict[str,Dict[str,Path]],Path]:\n chosen_root=None; best={}\n for r0 in _ixi_candidate_roots():\n root=_maybe_extract_ixi(r0)\n files=list(root.rglob(\"*.nii\"))+list(root.rglob(\"*.nii.gz\"))\n if not files: continue\n subs: Dict[str,Dict[str,Path]]={}\n for p in files:\n core=_strip_nii_ext(p.name)\n low=core.lower()\n kind=None; prefix=None\n for pref,k in ((\"mwp1\",\"GM\"),(\"mwp2\",\"WM\"),(\"wj\",\"JAC\"),(\"wm\",\"MRI\")):\n if low.startswith(pref): kind=k; prefix=pref; break\n if kind is None: continue\n suffix=core[len(prefix):].strip(\"_- .\")\n # Some releases place generic mwp1.nii/mwp2.nii/wj.nii/wm.nii in subject folders.\n if not suffix:\n rel=p.parent.relative_to(root) if p.parent!=root else Path(p.parent.name)\n pid=str(rel).replace(os.sep,\"__\")\n else:\n pid=suffix\n # Restrict mixed OASIS/Epilepsy/IXI packages to IXI when the identifier exposes source.\n combined=(pid+\" \"+str(p.parent)).lower()\n if (\"oasis\" in combined or \"epilep\" in combined) and \"ixi\" not in combined:\n continue\n subs.setdefault(pid,{})[kind]=p\n complete={pid:d for pid,d in subs.items() if all(k in d for k in (\"GM\",\"WM\",\"JAC\",\"MRI\"))}\n print(f\"[IXI DISCOVERY] root={root}: NIfTI={len(files)}, complete four-map subjects={len(complete)}\",flush=True)\n if len(complete)>len(best): best=complete; chosen_root=root\n if len(complete)>=8 and \"ixi\" in str(root).lower(): break\n if len(best)<8:\n raise FileNotFoundError(f\"Preprocessed IXI MRI requires >=8 subjects with GM/WM/Jacobian/MRI maps; found {len(best)}\")\n return best,chosen_root\n\n\ndef load_ixi() -> DatasetBundle:\n subs,root=_ixi_map_subjects()\n search_ids,val_ids,test_ids=deterministic_split(list(subs),IXI_MAX_SUBJECTS)\n maps=(\"GM\",\"WM\",\"JAC\",\"MRI\")\n def subject(pid:str)->np.ndarray:\n arr=[]\n for k in maps:\n arr.append(crop_or_pad_zyx(_load_nii_zyx(subs[pid][k]),MRI_CROP_ZYX))\n return np.stack(arr,axis=0).astype(np.float32) # map,z,y,x\n def stack(ids:List[str])->np.ndarray: return np.stack([subject(i) for i in ids],axis=0)\n raw_search=stack(search_ids); raw_val=stack(val_ids)\n mu,sd=channel_stats(raw_search,1,True)\n search=normalize_channels(raw_search,1,mu,sd); valid=normalize_channels(raw_val,1,mu,sd)\n del raw_search,raw_val; gc.collect()\n tests=[]\n for i in range(0,len(test_ids),max(1,IXI_GROUP_SIZE)):\n g=test_ids[i:i+max(1,IXI_GROUP_SIZE)]\n if not g: continue\n x=normalize_channels(stack(g),1,mu,sd)\n tests.append(TestUnit((\"subjects_\" if len(g)>1 else \"subject_\")+\"_\".join(g),x,None))\n progress_message(\"IXI SPLIT\",\"Subject-disjoint four-map MRI extension configured\",\n f\"development subjects={len(search_ids)}; validation={len(val_ids)}; held-out subjects={len(test_ids)} in {len(tests)} groups; maps={maps}; crop={MRI_CROP_ZYX}\")\n return DatasetBundle(\"ixi\",\"Preprocessed IXI MRI\",(\"P\",\"M\",\"Z\",\"Y\",\"X\"),1,maps,\n search,valid,tests,mu,sd,True,str(root),\"MRI-derived map intensity\")\n\n\ndef verify_all_five_dataset_sources() -> None:\n \"\"\"Fail before scientific measurements unless all five configured sources are discoverable.\"\"\"\n progress_message(\"DATASET GATE\", \"Verifying all five attached datasets before scientific work\",\n \"AMS/GEFS, MR-MS, INSTANT-ODC, HyperLeaf2024, IXI\")\n rows=[]; errors=[]\n def record(name, ok, detail):\n rows.append({\"dataset\":name,\"status\":\"OK\" if ok else \"FAIL\",\"detail\":str(detail)})\n print(f\"[DATASET GATE] {name}: {'OK' if ok else 'FAIL'} — {detail}\", flush=True)\n if not ok: errors.append(f\"{name}: {detail}\")\n try:\n tr=_find_ams_archive(\"train\"); te=_find_ams_archive(\"test\")\n record(\"ams_gefs\", True, f\"train={tr.name}; test={te.name}\")\n except Exception as e: record(\"ams_gefs\", False, e)\n try:\n pats=_find_mrms_patients(); record(\"mrms\", len(pats)>=8, f\"complete patients={len(pats)}\")\n except Exception as e: record(\"mrms\", False, e)\n try:\n subs=_instant_odc_subjects(); record(\"instant_odc\", len(subs)>=8, f\"complete subjects={len(subs)}\")\n except Exception as e: record(\"instant_odc\", False, e)\n try:\n root=_find_hyperleaf_root(); plots,_=_hyperleaf_plot_groups(root)\n record(\"hyperleaf\", len(plots)>=8, f\"root={root.name}; usable plots={len(plots)}\")\n except Exception as e: record(\"hyperleaf\", False, e)\n try:\n subs,root=_ixi_map_subjects(); record(\"ixi\", len(subs)>=8, f\"root={root.name if root else 'unknown'}; complete subjects={len(subs)}\")\n except Exception as e: record(\"ixi\", False, e)\n pd.DataFrame(rows).to_csv(OUT/\"dataset_validation_gate.csv\", index=False)\n if errors:\n raise RuntimeError(\"FIVE-DATASET GATE FAILED before scientific measurements: \" + \" | \".join(errors))\n progress_message(\"DATASET GATE\", \"ALL FIVE DATASETS VALID\", \"Scientific measurements may start.\")\n\n# -----------------------------------------------------------------------------\n# Compressor adapters\n# -----------------------------------------------------------------------------\n@dataclass\nclass CompressionResult:\n payload_bytes: int\n reconstruction: np.ndarray\n encode_seconds: float\n decode_seconds: float\n detail: str = \"\"\n\nclass Adapter:\n name = \"base\"\n max_order = 1\n execution_target = \"cpu\"\n def version(self) -> str: return \"unknown\"\n def compress(self, a: np.ndarray, eps: float, work: Path) -> CompressionResult: raise NotImplementedError\n\nclass SZ3Adapter(Adapter):\n name=\"SZ3\"; max_order=4\n def __init__(self):\n from pysz import sz, szConfig, szErrorBoundMode\n self.sz=sz; self.szConfig=szConfig; self.szErrorBoundMode=szErrorBoundMode\n def version(self): return package_version(\"pysz\")\n def compress(self,a,eps,work):\n a=np.ascontiguousarray(a,dtype=np.float32); cfg=self.szConfig(); cfg.errorBoundMode=self.szErrorBoundMode.ABS; cfg.absErrorBound=float(eps)\n t0=time.perf_counter(); comp,_=self.sz.compress(a,cfg); t1=time.perf_counter()\n rec,_=self.sz.decompress(comp,np.float32,a.shape); t2=time.perf_counter()\n n=len(comp) if isinstance(comp,(bytes,bytearray,memoryview)) else int(np.asarray(comp).nbytes)\n return CompressionResult(n,np.asarray(rec,dtype=np.float32).reshape(a.shape),t1-t0,t2-t1)\n\n\nclass HPEZAdapter(Adapter):\n name=\"HPEZ\"; max_order=4\n def __init__(self,binpath:Path,commit:str): self.bin=binpath; self.commit=commit\n def version(self): return self.commit or \"unknown\"\n def compress(self,a,eps,work):\n a=np.ascontiguousarray(a,dtype=np.float32); sub=work/(\"hpez_\"+hashlib.md5((str(a.shape)+str(eps)).encode()).hexdigest()); shutil.rmtree(sub,ignore_errors=True); sub.mkdir(parents=True)\n raw=sub/\"input.f32\"; a.tofile(raw); dimflag=f\"-{a.ndim}\"; dimargs=[str(int(x)) for x in reversed(a.shape)]\n before=set(sub.iterdir()); cmd=[str(self.bin),\"-z\",\"-f\",\"-i\",str(raw),dimflag]+dimargs+[\"-M\",\"ABS\",\"-A\",str(float(eps)),\"-q\",\"4\"]\n t0=time.perf_counter(); p=sh(cmd,cwd=sub,timeout=1200); t1=time.perf_counter()\n if p.returncode!=0: raise RuntimeError(p.stdout[-2000:])\n created=[q for q in sub.iterdir() if q not in before and q!=raw and q.is_file() and q.stat().st_size None:\n pkgs=[\"pysz==1.0.3\",\"nibabel\",\"netCDF4\",\"psutil\",\"scikit-image\",\"matplotlib\",\"tifffile\",\"cmake<4\"]\n pinned={\"pysz\":\"1.0.3\"}\n required=(\"pysz\",\"nibabel\",\"netCDF4\",\"psutil\",\"scikit-image\",\"matplotlib\",\"tifffile\",\"cmake\")\n versions={name:package_version(name) for name in required}\n bad=[name for name in required if versions.get(name) in (\"unavailable\", \"unknown\", \"\")]\n bad += [name for name,want in pinned.items() if versions.get(name) != want]\n if not bad and not FORCE_REBUILD:\n progress_message(\"1/4 PYTHON PACKAGES\", \"VERIFIED — installation skipped\",\n \"; \".join(f\"{k}={v}\" for k,v in versions.items()))\n emit_event(\"software_verify\",\"SKIP\",detail=\"Python dependencies already satisfy pinned requirements\")\n return\n progress_message(\"1/4 PYTHON PACKAGES\", \"Installing/repairing Python dependencies\",\n \"Verification failed for: \" + \", \".join(sorted(set(bad))) + \". Packages: \" + \", \".join(pkgs))\n p=run_command_progress([sys.executable,\"-m\",\"pip\",\"install\",\"--progress-bar\",\"off\"]+pkgs,\n \"Python dependency installation\", timeout=1800, heartbeat=12)\n if p.returncode!=0:\n raise RuntimeError(p.stdout[-4000:])\n versions = {name: package_version(name) for name in required}\n progress_message(\"1/4 PYTHON PACKAGES\", \"Dependency installation complete\",\n \"; \".join(f\"{k}={v}\" for k,v in versions.items()))\n\n\ndef build_hpez() -> Tuple[Optional[Path],str,str]:\n src=TOOLS/\"QoZ\"; inst=TOOLS/\"hpez_install\"\n try:\n progress_message(\"3/4 HPEZ/QoZ\", \"Verifying existing HPEZ/QoZ before any build\", f\"Pinned commit: {HPEZ_COMMIT}; build jobs={BUILD_JOBS}\")\n if src.exists():\n current=sh([\"git\",\"rev-parse\",\"HEAD\"],cwd=src,timeout=20).stdout.strip()\n cands=list(inst.rglob(\"qoz\"))+list(inst.rglob(\"hpez\")) if inst.exists() else []\n binp=next((q for q in cands if q.is_file() and os.access(q,os.X_OK)),None)\n if current == HPEZ_COMMIT and binp is not None and not FORCE_REBUILD:\n progress_message(\"3/4 HPEZ/QoZ\", \"VERIFIED — build skipped\", f\"commit={current[:12]}; executable={binp}\")\n emit_event(\"software_verify\",\"SKIP\",compressor=\"HPEZ\",detail=f\"verified pinned executable {binp}\")\n return binp,current,\"verified existing build\"\n progress_message(\"3/4 HPEZ/QoZ\", \"Existing build is incomplete/stale\", f\"commit={current[:12] if current else 'unknown'}; executable={binp}; rebuilding required\")\n if not src.exists():\n run_command_progress([\"git\",\"clone\",HPEZ_REPO,str(src)], \"HPEZ clone\", timeout=300, check=True)\n else:\n progress_message(\"3/4 HPEZ/QoZ\", \"Source present; clone skipped\", str(src))\n run_command_progress([\"git\",\"fetch\",\"--all\"], \"HPEZ git fetch\", cwd=src, timeout=180, check=True)\n run_command_progress([\"git\",\"checkout\",HPEZ_COMMIT], \"HPEZ checkout pinned commit\", cwd=src, timeout=60, check=True)\n commit=sh([\"git\",\"rev-parse\",\"HEAD\"],cwd=src,check=True).stdout.strip()\n build=src/\"build\"; build.mkdir(exist_ok=True)\n run_command_progress([\"cmake\",f\"-DCMAKE_INSTALL_PREFIX={inst}\",\"..\"], \"HPEZ CMake configure\", cwd=build, timeout=600, check=True)\n run_command_progress([\"cmake\",\"--build\",\".\",\"-j\",str(BUILD_JOBS)], \"HPEZ compile\", cwd=build, timeout=3600, heartbeat=15, check=True)\n run_command_progress([\"cmake\",\"--install\",\".\"], \"HPEZ install\", cwd=build, timeout=600, check=True)\n cands=list(inst.rglob(\"qoz\"))+list(inst.rglob(\"hpez\"))\n binp=next((q for q in cands if q.is_file() and os.access(q,os.X_OK)),None)\n progress_message(\"3/4 HPEZ/QoZ\", \"HPEZ/QoZ setup complete\" if binp else \"HPEZ/QoZ binary not found\", f\"commit={commit[:12]}; binary={binp}\")\n return binp,commit,\"ok\" if binp else \"binary not found\"\n except Exception as e:\n progress_message(\"3/4 HPEZ/QoZ\", \"HPEZ/QoZ setup failed\", str(e))\n return None,\"\",str(e)\n\n\n\n\n\n\n\n\n\n\n\ndef setup_adapters() -> Tuple[List[Adapter],pd.DataFrame]:\n progress_message(\"SOFTWARE SETUP\",\"Preparing two-codec CPU benchmark\",\n \"Required: SZ3 and HPEZ. No GPU, SPERR, SZx, vecSZ, or MGARD build is attempted.\")\n emit_event(\"software_setup\",\"START\",detail=\"SZ3,HPEZ\")\n pip_install(); rows=[]; adapters=[]\n try:\n a=SZ3Adapter(); adapters.append(a)\n rows.append({\"compressor\":a.name,\"available\":1,\"version\":a.version(),\"execution_target\":\"cpu\",\"detail\":\"verified Python package\"})\n except Exception as e:\n rows.append({\"compressor\":\"SZ3\",\"available\":0,\"version\":\"\",\"execution_target\":\"cpu\",\"detail\":str(e)})\n hb,hcommit,hstatus=build_hpez()\n rows.append({\"compressor\":\"HPEZ\",\"available\":int(hb is not None),\"version\":hcommit,\"execution_target\":\"cpu\",\"detail\":hstatus})\n if hb is not None: adapters.append(HPEZAdapter(hb,hcommit))\n # Keep output order stable for the paper.\n order={n:i for i,n in enumerate(FINAL_CODECS)}\n adapters.sort(key=lambda a: order.get(a.name,999))\n df=pd.DataFrame(rows); df.to_csv(OUT/\"software_status.csv\",index=False)\n for _,r in df.iterrows():\n emit_event(\"software_setup\",\"OK\" if int(r.available)==1 else \"FAIL\",compressor=str(r.compressor),detail=f\"target=cpu; {r.detail}\")\n missing=[n for n in REQUIRED_CODECS if not ((df.compressor==n)&(df.available==1)).any()]\n progress_message(\"SOFTWARE CHECK\",\"CPU codec setup summary\",\", \".join(f\"{r.compressor}={'OK' if int(r.available) else 'FAIL'}[cpu]\" for _,r in df.iterrows()))\n if STRICT_CODECS and missing: raise RuntimeError(\"required CPU compressors unavailable: \"+\", \".join(missing))\n return adapters,df\n\n# -----------------------------------------------------------------------------\n# Measurement and layout selection\n# -----------------------------------------------------------------------------\ndef _stream_error_stats(ref: np.ndarray, rec: np.ndarray, mask: Optional[np.ndarray] = None,\n buffer_values: int = 262144) -> Tuple[float, float, int]:\n \"\"\"Return sum of squared error, max abs error, and count with bounded temporaries.\"\"\"\n if ref.shape != rec.shape:\n raise ValueError(f\"shape mismatch: {ref.shape} vs {rec.shape}\")\n ops = [ref, rec]\n op_flags = [[\"readonly\"], [\"readonly\"]]\n op_dtypes = [np.float32, np.float32]\n if mask is not None:\n mask = np.broadcast_to(np.asarray(mask, dtype=bool), ref.shape)\n ops.append(mask); op_flags.append([\"readonly\"]); op_dtypes.append(np.bool_)\n it = np.nditer(ops, flags=[\"external_loop\", \"buffered\", \"zerosize_ok\"],\n op_flags=op_flags, op_dtypes=op_dtypes, buffersize=buffer_values, order=\"K\")\n ss = 0.0; ma = 0.0; count = 0\n for vals in it:\n a, b = vals[0], vals[1]\n d = np.asarray(a, dtype=np.float32) - np.asarray(b, dtype=np.float32)\n if mask is not None:\n m = np.asarray(vals[2], dtype=bool)\n if not np.any(m):\n continue\n d = d[m]\n if d.size:\n d64 = d.astype(np.float64, copy=False)\n ss += float(np.dot(d64, d64))\n ma = max(ma, float(np.max(np.abs(d))))\n count += int(d.size)\n return ss, ma, count\n\ndef measure(ref:np.ndarray,rec:np.ndarray,dr:float)->Tuple[float,float,float,float,float]:\n ss, ma, n = _stream_error_stats(ref, rec)\n mse=max(ss/max(n,1),1e-30); rmse=math.sqrt(mse); nrmse=rmse/dr; psnr=20*math.log10(dr/rmse)\n return mse,rmse,nrmse,psnr,ma\n\n\ndef one_point(adapter:Adapter,x:np.ndarray,layout:Layout,eps:float,dr:float,work:Path)->Tuple[dict,np.ndarray]:\n a=fold_tensor(x,layout); r=adapter.compress(a,eps,work); rec=unfold_tensor(r.reconstruction,x.shape,layout)\n mse,rmse,nrmse,psnr,ma=measure(x,rec,dr); total=r.payload_bytes+COMMON_WRAPPER_BYTES\n enc_mbps = x.nbytes/1e6/r.encode_seconds if np.isfinite(r.encode_seconds) and r.encode_seconds>0 else np.nan\n dec_mbps = x.nbytes/1e6/r.decode_seconds if np.isfinite(r.decode_seconds) and r.decode_seconds>0 else np.nan\n row={\"payload_bytes\":r.payload_bytes,\"wrapper_bytes\":COMMON_WRAPPER_BYTES,\"total_bytes\":total,\n \"bits_per_value\":8.0*total/x.size,\"compression_ratio\":x.nbytes/max(total,1),\"mse\":mse,\"rmse\":rmse,\"nrmse\":nrmse,\"psnr_db\":psnr,\n \"max_abs_error\":ma,\"bound_ratio\":ma/max(eps,1e-30),\"encode_seconds\":r.encode_seconds,\"decode_seconds\":r.decode_seconds,\n \"encode_MBps\":enc_mbps,\"decode_MBps\":dec_mbps,\"execution_target\":adapter.execution_target,\"adapter_detail\":r.detail}\n return row,rec\n\n\ndef codec_fullshape_preflight(bundle:DatasetBundle,adapter:Adapter,dr:float)->Tuple[bool,str]:\n \"\"\"Exercise one real validation-shape encode/decode before layout search.\n\n Required CPU codecs fail the dataset if this test fails.\n The preflight result is diagnostic only and never selects a representation.\n \"\"\"\n m=min(adapter.max_order,bundle.d)\n lay=reference_layout(bundle.d,m)\n eps=eps_from_nominal(float(VALID_NOMINAL_GRID[len(VALID_NOMINAL_GRID)//2]),dr)\n shape=tuple(int(x) for x in fold_tensor(bundle.validation,lay).shape)\n detail=(f\"native_shape={tuple(bundle.validation.shape)}; presented_shape={shape}; \"\n f\"d={bundle.d}->m={m}; tolerance={eps:.8g}\")\n progress_message(\"PREFLIGHT\", f\"{bundle.dataset_id}/{adapter.name}: full-shape codec check\", detail)\n emit_event(\"codec_preflight\",\"START\",dataset=bundle.dataset_id,compressor=adapter.name,detail=detail)\n work=ROOT/\"tmp\"/bundle.dataset_id/(adapter.name.replace(\"/\",\"_\")+\"_preflight\")\n work.mkdir(parents=True,exist_ok=True)\n try:\n row,_=one_point(adapter,bundle.validation,lay,eps,dr,work)\n ratio=float(row.get(\"bound_ratio\", np.inf))\n if (not np.isfinite(ratio)) or ratio > BOUND_ACCEPT_RATIO:\n raise RuntimeError(f\"measured bound ratio {ratio} exceeds acceptance ratio {BOUND_ACCEPT_RATIO}\")\n except Exception as e:\n msg=(f\"{adapter.name} full-shape preflight failed before layout search. {detail}. \"\n f\"Backend error: {type(e).__name__}: {e}\")\n emit_event(\"codec_preflight\",\"FAIL\",dataset=bundle.dataset_id,compressor=adapter.name,detail=msg[:4000])\n if adapter.name in OPTIONAL_CODECS:\n progress_message(\n \"CODEC EXCLUDED\",\n f\"{bundle.dataset_id}/{adapter.name}: excluded from this dataset\",\n \"The codec could not verify the common fixed point-wise error contract on the full folded tensor. \"\n \"It will not enter layout selection, RD summaries, external envelopes, or paper aggregates for this dataset.\\n\"\n f\"Reason: {msg}\"\n )\n return False,msg\n raise RuntimeError(msg) from e\n okmsg=(f\"{detail}; bpv={row['bits_per_value']:.6g}; measured_PSNR={row['psnr_db']:.3f}; \"\n f\"bound_ratio={row['bound_ratio']:.8f}\")\n emit_event(\"codec_preflight\",\"OK\",dataset=bundle.dataset_id,compressor=adapter.name,detail=okmsg)\n progress_message(\"PREFLIGHT OK\",f\"{bundle.dataset_id}/{adapter.name}: common error contract verified\",\n f\"bpv={row['bits_per_value']:.6g}; measured PSNR={row['psnr_db']:.3f} dB; bound ratio={row['bound_ratio']:.8f}\")\n del row\n gc.collect()\n return True,okmsg\n\n\ndef _pareto_rd_points(rows: Sequence[dict]) -> pd.DataFrame:\n d=pd.DataFrame(list(rows))\n if len(d)<2: return pd.DataFrame()\n d=d.copy(); d[\"psnr_db\"]=pd.to_numeric(d[\"psnr_db\"],errors=\"coerce\"); d[\"bits_per_value\"]=pd.to_numeric(d[\"bits_per_value\"],errors=\"coerce\")\n d=d[np.isfinite(d.psnr_db)&np.isfinite(d.bits_per_value)].sort_values([\"psnr_db\",\"bits_per_value\"]).reset_index(drop=True)\n if len(d)<2: return d\n keep=[]\n for i,r in d.iterrows():\n dominated=((d.psnr_db >= float(r.psnr_db)-1e-12) & (d.bits_per_value <= float(r.bits_per_value)+1e-12) &\n ((d.psnr_db > float(r.psnr_db)+1e-12) | (d.bits_per_value < float(r.bits_per_value)-1e-12))).any()\n if not dominated: keep.append(i)\n return d.loc[keep].sort_values(\"psnr_db\").reset_index(drop=True)\n\ndef _interp_layout_rate(d: pd.DataFrame, psnr: float) -> Optional[float]:\n q=_pareto_rd_points(d.to_dict(\"records\")) if isinstance(d,pd.DataFrame) else _pareto_rd_points(d)\n if len(q)<2: return None\n x=q.psnr_db.to_numpy(float); y=q.bits_per_value.to_numpy(float)\n if psnr < float(x.min())-1e-9 or psnr > float(x.max())+1e-9: return None\n return float(np.interp(psnr,x,y))\n\ndef validation_layouts(bundle:DatasetBundle,adapter:Adapter,dr:float,all_search_rows:List[dict])->Tuple[Layout,List[Layout]]:\n m=min(adapter.max_order,bundle.d); layouts=enumerate_layouts(bundle.d,m); ref=reference_layout(bundle.d,m)\n log(f\"[{bundle.dataset_id}/{adapter.name}] layout space d={bundle.d}->m={m}: {len(layouts)}\")\n emit_event(\"layout_search\",\"START\",dataset=bundle.dataset_id,compressor=adapter.name,detail=f\"{bundle.d}->{m}; {len(layouts)} layouts; selection_revision={SELECTION_REVISION}\")\n # Stage 1: evaluate the locality proxy over the complete permutation/folding space.\n # Scores are shared across codecs that expose the same presented tensor order.\n proxy_key=(bundle.dataset_id,m)\n if proxy_key in PROXY_LAYOUT_CACHE:\n proxy_rows=PROXY_LAYOUT_CACHE[proxy_key]\n log(f\"[PROXY REUSE] {bundle.dataset_id}/{adapter.name}: reusing {len(proxy_rows)} codec-independent scores for presented order m={m}\")\n else:\n proxy_rows = proxy_layout_scores_fast(\n bundle.search, layouts, bundle.axis_names, bundle.dataset_id, f\"shared-m{m}\"\n )\n PROXY_LAYOUT_CACHE[proxy_key]=proxy_rows\n for sc, lay in proxy_rows:\n all_search_rows.append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"stage\":\"proxy\",\"layout\":lay.machine_label(),\"layout_semantic\":lay.label(bundle.axis_names),\"score\":sc})\n # Uniform shortlist: retain strong global proxy layouts and strong layouts within\n # every contiguous folding partition. This avoids dataset-specific candidate injection\n # while preventing a single global proxy ranking from excluding an entire folding pattern.\n candidate_map: Dict[str, Layout] = {}\n for _, lay in proxy_rows[:PROXY_GLOBAL_KEEP]:\n candidate_map[lay.machine_label()] = lay\n by_partition: Dict[Tuple[int, ...], List[Tuple[float, Layout]]] = {}\n for sc, lay in proxy_rows:\n by_partition.setdefault(tuple(lay.cuts), []).append((sc, lay))\n for cuts in sorted(by_partition):\n ranked_partition = sorted(by_partition[cuts], key=lambda z: (z[0], z[1].machine_label()))\n for _, lay in ranked_partition[:PROXY_PER_PARTITION_KEEP]:\n candidate_map[lay.machine_label()] = lay\n candidates = [candidate_map[k] for k in sorted(candidate_map)]\n log(f\"[PROXY FINALISTS] {bundle.dataset_id}/{adapter.name}: retaining {len(candidates)} of {len(layouts)} layouts \"\n f\"(global top={PROXY_GLOBAL_KEEP}; per-fold-partition top={PROXY_PER_PARTITION_KEEP}; partitions={len(by_partition)})\")\n for rank_i, (sc, lay) in enumerate(proxy_rows[:min(10, len(proxy_rows))], 1):\n log(f\" global proxy rank {rank_i:2d}: {lay.label(bundle.axis_names):20s} score={sc:.6f}\")\n if ref not in candidates:\n candidates.append(ref)\n probe=center_crop_probe(bundle.search); work=ROOT/\"tmp\"/bundle.dataset_id/adapter.name.replace(\"/\",\"_\"); work.mkdir(parents=True,exist_ok=True)\n coarse=[]; per_quality=[]\n log(f\"[COARSE SEARCH] {bundle.dataset_id}/{adapter.name}: \"\n f\"{len(candidates)} proxy finalists x {len(SEARCH_NOMINAL_GRID)} tolerances\")\n coarse_t0=time.perf_counter()\n for ci, lay in enumerate(candidates, 1):\n rates=[]; ok=True\n log(f\"[COARSE {ci}/{len(candidates)}] {bundle.dataset_id}/{adapter.name}: \"\n f\"layout={lay.label(bundle.axis_names)}\")\n for db in SEARCH_NOMINAL_GRID:\n eps=eps_from_nominal(db,dr)\n try:\n row,_=one_point(adapter,probe,lay,eps,dr,work); rates.append(row[\"bits_per_value\"]); per_quality.append((float(db),float(row[\"bits_per_value\"]),lay))\n except Exception as e:\n ok=False; all_search_rows.append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"stage\":\"coarse_error\",\"layout\":lay.machine_label(),\"layout_semantic\":lay.label(bundle.axis_names),\"score\":np.inf,\"error\":str(e)[:500]}); break\n if ok:\n score=float(np.mean(np.log(np.maximum(rates,1e-12)))); coarse.append((score,lay)); all_search_rows.append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"stage\":\"coarse\",\"layout\":lay.machine_label(),\"layout_semantic\":lay.label(bundle.axis_names),\"score\":score})\n elapsed=time.perf_counter()-coarse_t0\n best_now=min(coarse, key=lambda z:(z[0], z[1].machine_label()))\n log(f\" coarse score={score:.6f}; rates=\" + \", \".join(f\"{r:.4f}\" for r in rates))\n log(f\" current coarse best={best_now[1].label(bundle.axis_names)} \"\n f\"score={best_now[0]:.6f}; elapsed={elapsed:.1f}s\")\n if not coarse: raise RuntimeError(f\"no coarse layout succeeded for {adapter.name}\")\n coarse.sort(key=lambda z:(z[0], z[1].machine_label())); finalists=[z[1] for z in coarse[:COARSE_KEEP]]\n # Retain the strongest layouts at each coarse tolerance in addition to the aggregate ranking.\n for db in SEARCH_NOMINAL_GRID:\n q=sorted(((r, lay.machine_label(), lay) for d,r,lay in per_quality if abs(d-float(db))<1e-9), key=lambda z:(z[0], z[1]))\n for _,_,lay in q[:2]:\n if lay not in finalists: finalists.append(lay)\n if ref not in finalists:\n finalists.append(ref)\n validation_errors=[]; curves: Dict[str,List[dict]]={}\n log(f\"[VALIDATION SEARCH] {bundle.dataset_id}/{adapter.name}: \"\n f\"{len(finalists)} finalists x {len(VALID_NOMINAL_GRID)} full-validation tolerances\")\n val_t0=time.perf_counter()\n for vi, lay in enumerate(finalists, 1):\n rows=[]; ok=True\n log(f\"[VALIDATION {vi}/{len(finalists)}] {bundle.dataset_id}/{adapter.name}: \"\n f\"layout={lay.label(bundle.axis_names)}\")\n for db in VALID_NOMINAL_GRID:\n eps=eps_from_nominal(db,dr)\n try:\n row,_=one_point(adapter,bundle.validation,lay,eps,dr,work)\n rows.append({\"nominal_psnr\":float(db),\"bits_per_value\":row[\"bits_per_value\"],\"psnr_db\":row[\"psnr_db\"],\"bound_ratio\":row[\"bound_ratio\"]})\n log(f\" nominal={float(db):g}: bpv={row['bits_per_value']:.6f}; \"\n f\"measured_PSNR={row['psnr_db']:.3f} dB; bound={row['bound_ratio']:.8f}\")\n all_search_rows.append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"stage\":\"validation_point\",\"layout\":lay.machine_label(),\"layout_semantic\":lay.label(bundle.axis_names),\"nominal_psnr\":float(db),\"bits_per_value\":row[\"bits_per_value\"],\"psnr_db\":row[\"psnr_db\"],\"bound_ratio\":row[\"bound_ratio\"]})\n except Exception as e:\n ok=False; errmsg=f\"{type(e).__name__}: {e}\"; validation_errors.append((lay.label(bundle.axis_names),float(db),errmsg))\n all_search_rows.append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"stage\":\"validation_error\",\"layout\":lay.machine_label(),\"layout_semantic\":lay.label(bundle.axis_names),\"score\":np.inf,\"error\":errmsg[:1000]})\n log(f\"VALIDATION ERROR {bundle.dataset_id}/{adapter.name} layout={lay.label(bundle.axis_names)} nominal={db:g}: {errmsg}\"); break\n if ok and len(rows)>=2:\n curves[lay.machine_label()]=rows\n elapsed=time.perf_counter()-val_t0\n log(f\" validation curve complete ({vi}/{len(finalists)}); elapsed={elapsed:.1f}s\")\n if not curves:\n samples=\" | \".join(f\"layout={a}, nominal={b:g}: {c}\" for a,b,c in validation_errors[:3]) or \"no exception detail captured\"\n raise RuntimeError(f\"no validation layout succeeded for {adapter.name}. Representative errors: {samples}\")\n clean={k:_pareto_rd_points(v) for k,v in curves.items()}\n clean={k:v for k,v in clean.items() if len(v)>=2}\n lo=max(float(v.psnr_db.min()) for v in clean.values()); hi=min(float(v.psnr_db.max()) for v in clean.values())\n if hi<=lo: raise RuntimeError(f\"validation layouts for {adapter.name} have no common measured-PSNR overlap\")\n counts={k:0 for k in clean}; regrets={k:[] for k in clean}\n for p in np.linspace(lo,hi,VALID_WINNER_GRID):\n opts=[]\n for k,v in clean.items():\n r=_interp_layout_rate(v,p)\n if r is not None: opts.append((r,k))\n if not opts: continue\n opts.sort(); best_r,best_k=opts[0]; counts[best_k]+=1\n for r,k in opts: regrets[k].append(100.0*(r-best_r)/max(best_r,1e-30))\n rank=[]\n for k in clean:\n mr=float(np.mean(regrets[k])) if regrets[k] else np.inf\n rank.append((-int(counts[k]),mr,k))\n lay=parse_layout_machine_label(k)\n all_search_rows.append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"stage\":\"validation_rank\",\"layout\":k,\"layout_semantic\":lay.label(bundle.axis_names),\"minimum_rate_grid_count\":int(counts[k]),\"mean_relative_excess_rate_pct\":mr,\"common_psnr_lo\":lo,\"common_psnr_hi\":hi})\n rank.sort(); ranked=[parse_layout_machine_label(k) for _,_,k in rank]\n selected=ranked[0]\n emit_event(\"layout_search\",\"OK\",dataset=bundle.dataset_id,compressor=adapter.name,detail=f\"selected {selected.label(bundle.axis_names)} from {len(ranked)} validated finalists using the uniform measured-PSNR criterion\")\n progress_message(\"LAYOUT SELECT\",f\"{bundle.dataset_id}/{adapter.name}: validation-selected layout\",f\"selected={selected.label(bundle.axis_names)}; minimum-rate grid count={counts[selected.machine_label()]}/{VALID_WINNER_GRID}; mean relative excess rate={float(np.mean(regrets[selected.machine_label()])):.4f}%; validation PSNR overlap={lo:.2f}..{hi:.2f} dB\")\n return selected,ranked\n\n\ndef prune_stale_selected_measurements(tables: Dict[str,List[dict]], dataset: str, compressor: str, selected_layout: str) -> None:\n \"\"\"Keep reusable exact-layout points while removing selected-policy rows from an older selection revision.\"\"\"\n for name in (\"raw_rd_points\",\"per_channel_metrics\",\"geometry_metrics\"):\n kept=[]\n for r in tables.get(name,[]):\n if str(r.get(\"dataset\"))==dataset and str(r.get(\"compressor\"))==compressor and str(r.get(\"layout_kind\"))==\"selected\":\n if name==\"geometry_metrics\":\n continue # geometry rows do not store the layout label; recompute cheaply\n if str(r.get(\"layout\", \"\")) != selected_layout:\n continue\n kept.append(r)\n tables[name]=kept\n\ndef geometry_rows(bundle:DatasetBundle,unit:TestUnit,layout:Layout,kind:str,compressor:str)->List[dict]:\n # Geometry diagnostics are descriptive, so estimate them on a deterministic\n # bounded center crop rather than duplicating a full held-out tensor.\n probe = center_crop_probe(unit.x, GEOMETRY_MAX_VALUES)\n y=fold_tensor(probe,layout); scale=float(np.std(y,dtype=np.float64))+1e-12; rows=[]\n for ax,n in enumerate(y.shape):\n val=float(np.mean(np.abs(np.diff(y,axis=ax)),dtype=np.float64))/scale if n>1 else np.nan\n rows.append({\"dataset\":bundle.dataset_id,\"unit\":unit.unit_id,\"compressor\":compressor,\"layout_kind\":kind,\"presented_axis\":ax,\"presented_size\":n,\"geometry_sample_values\":int(probe.size),\"normalized_mean_abs_first_difference\":val})\n return rows\n\n\ndef channel_metric_rows(bundle:DatasetBundle,unit:TestUnit,rec:np.ndarray,base:dict,layout_kind:str)->List[dict]:\n rows=[]; ca=bundle.channel_axis\n nchan=len(bundle.channel_names)\n if nchan <= MAX_DIAGNOSTIC_CHANNELS:\n channel_indices=list(range(nchan))\n else:\n channel_indices=sorted(set(int(round(x)) for x in np.linspace(0,nchan-1,MAX_DIAGNOSTIC_CHANNELS)))\n for c in channel_indices:\n name=bundle.channel_names[c]\n sl=[slice(None)]*unit.x.ndim; sl[ca]=c\n refc=unit.x[tuple(sl)]; recc=rec[tuple(sl)]\n ss,ma,n=_stream_error_stats(refc,recc); rmse=math.sqrt(ss/max(n,1))\n native_rmse=rmse*float(bundle.normalization_std[c]); native_ma=ma*float(bundle.normalization_std[c])\n rr={k:base[k] for k in (\"dataset\",\"unit\",\"compressor\",\"nominal_psnr\",\"absolute_tolerance\")}; rr.update({\"layout_kind\":layout_kind,\"channel\":name,\"rmse_normalized\":rmse,\"max_abs_error_normalized\":ma,\"rmse_native\":native_rmse,\"max_abs_error_native\":native_ma,\"native_unit_label\":bundle.native_unit_label})\n # ROI mask is P,Z,Y,X for MRI and broadcasts over visit/modality axes.\n if unit.roi_mask is not None:\n try:\n mask=np.asarray(unit.roi_mask,dtype=bool)\n while mask.ndim < refc.ndim: mask=np.expand_dims(mask,axis=1)\n mask=np.broadcast_to(mask,refc.shape)\n rss,_,rn=_stream_error_stats(refc,recc,mask)\n rr[\"roi_rmse_normalized\"]=math.sqrt(rss/rn) if rn else np.nan\n rr[\"roi_rmse_native\"]=rr[\"roi_rmse_normalized\"]*float(bundle.normalization_std[c]) if np.isfinite(rr[\"roi_rmse_normalized\"]) else np.nan\n except Exception: rr[\"roi_rmse_normalized\"]=np.nan; rr[\"roi_rmse_native\"]=np.nan\n rows.append(rr)\n return rows\n\n\ndef integrated_saving(curve_opt:pd.DataFrame,curve_ref:pd.DataFrame,n: int=101)->Optional[float]:\n a=curve_opt.sort_values(\"psnr_db\"); b=curve_ref.sort_values(\"psnr_db\")\n if len(a)<2 or len(b)<2: return None\n lo=max(float(a.psnr_db.min()),float(b.psnr_db.min())); hi=min(float(a.psnr_db.max()),float(b.psnr_db.max()))\n if hi<=lo: return None\n p=np.linspace(lo,hi,n); ra=np.interp(p,a.psnr_db,a.bits_per_value); rb=np.interp(p,b.psnr_db,b.bits_per_value)\n delta=float(np.mean(np.log(np.maximum(ra,1e-30))-np.log(np.maximum(rb,1e-30))))\n return 100.0*(1.0-math.exp(delta))\n\n\ndef bootstrap_ci(vals:Sequence[float],nboot:int=4000)->Tuple[float,float]:\n x=np.asarray(vals,dtype=float); x=x[np.isfinite(x)]\n if len(x)==0:return np.nan,np.nan\n rng=np.random.default_rng(SEED+7); means=np.empty(nboot)\n for i in range(nboot): means[i]=np.mean(rng.choice(x,size=len(x),replace=True))\n return float(np.percentile(means,2.5)),float(np.percentile(means,97.5))\n\n\n\ndef _measured_psnr_intervals() -> List[Tuple[float, float]]:\n return [\n (PSNR_INTERVAL_START_DB + i * PSNR_INTERVAL_WIDTH_DB,\n PSNR_INTERVAL_START_DB + (i + 1) * PSNR_INTERVAL_WIDTH_DB)\n for i in range(PSNR_INTERVAL_COUNT)\n ]\n\n\ndef _integrated_saving_interval(curve_opt: pd.DataFrame, curve_ref: pd.DataFrame,\n band_lo: float, band_hi: float) -> Optional[dict]:\n a = _pareto_rd_points(curve_opt.to_dict(\"records\"))\n b = _pareto_rd_points(curve_ref.to_dict(\"records\"))\n if len(a) < 2 or len(b) < 2:\n return None\n support_lo = max(float(a.psnr_db.min()), float(b.psnr_db.min()))\n support_hi = min(float(a.psnr_db.max()), float(b.psnr_db.max()))\n lo = max(float(band_lo), support_lo)\n hi = min(float(band_hi), support_hi)\n coverage = hi - lo\n if coverage < PSNR_INTERVAL_MIN_COVERAGE_DB:\n return None\n p = np.linspace(lo, hi, PSNR_INTERVAL_GRID_POINTS)\n rs = np.interp(p, a.psnr_db.to_numpy(float), a.bits_per_value.to_numpy(float))\n rr = np.interp(p, b.psnr_db.to_numpy(float), b.bits_per_value.to_numpy(float))\n delta = float(np.mean(np.log(np.maximum(rs, 1e-30)) - np.log(np.maximum(rr, 1e-30))))\n return {\n \"saving_pct\": 100.0 * (1.0 - math.exp(delta)),\n \"effective_psnr_lo\": lo,\n \"effective_psnr_hi\": hi,\n \"coverage_db\": coverage,\n \"curve_support_lo\": support_lo,\n \"curve_support_hi\": support_hi,\n }\n\n\ndef measured_psnr_interval_tables(valid: pd.DataFrame) -> Tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]:\n \"\"\"Compute the same fixed measured-PSNR interval analysis for every dataset/codec.\n\n No interval participates in layout selection. The calculation is a descriptive\n decomposition of held-out reference and selected RD curves after all layouts are fixed.\n \"\"\"\n per_unit = []\n intervals = _measured_psnr_intervals()\n for (ds, unit, codec), g in valid.groupby([\"dataset\", \"unit\", \"compressor\"]):\n sel = g[g.layout_kind.astype(str).eq(\"selected\")]\n ref = g[g.layout_kind.astype(str).eq(\"reference\")]\n for band_index, (blo, bhi) in enumerate(intervals):\n r = _integrated_saving_interval(sel, ref, blo, bhi)\n if r is None:\n continue\n per_unit.append({\n \"dataset\": ds, \"unit\": unit, \"compressor\": codec,\n \"interval_index\": band_index,\n \"interval_lo_db\": blo, \"interval_hi_db\": bhi,\n \"interval_label\": f\"{blo:g}-{bhi:g} dB\",\n \"rate_reduction_pct\": r[\"saving_pct\"],\n \"effective_psnr_lo\": r[\"effective_psnr_lo\"],\n \"effective_psnr_hi\": r[\"effective_psnr_hi\"],\n \"coverage_db\": r[\"coverage_db\"],\n \"curve_support_lo\": r[\"curve_support_lo\"],\n \"curve_support_hi\": r[\"curve_support_hi\"],\n })\n pdf = pd.DataFrame(per_unit)\n summary_rows = []\n if len(pdf):\n for (ds, codec, idx, blo, bhi, label), g in pdf.groupby(\n [\"dataset\", \"compressor\", \"interval_index\", \"interval_lo_db\", \"interval_hi_db\", \"interval_label\"],\n sort=True,\n ):\n x = g.rate_reduction_pct.to_numpy(float)\n ci_lo, ci_hi = bootstrap_ci(x)\n summary_rows.append({\n \"dataset\": ds, \"compressor\": codec, \"interval_index\": int(idx),\n \"interval_lo_db\": float(blo), \"interval_hi_db\": float(bhi), \"interval_label\": label,\n \"test_units\": int(len(x)), \"mean_rate_reduction_pct\": float(np.mean(x)),\n \"median_rate_reduction_pct\": float(np.median(x)),\n \"min_rate_reduction_pct\": float(np.min(x)), \"max_rate_reduction_pct\": float(np.max(x)),\n \"std_rate_reduction_pct\": float(np.std(x, ddof=1)) if len(x) > 1 else 0.0,\n \"bootstrap95_lo\": ci_lo, \"bootstrap95_hi\": ci_hi,\n \"positive_units\": int(np.sum(x > 0)),\n \"minimum_coverage_db\": float(g.coverage_db.min()),\n \"mean_coverage_db\": float(g.coverage_db.mean()),\n })\n sdf = pd.DataFrame(summary_rows)\n\n shape_rows = []\n if len(sdf):\n for (ds, codec), g in sdf.groupby([\"dataset\", \"compressor\"], sort=True):\n q = g.sort_values(\"interval_index\").reset_index(drop=True)\n vals = q.mean_rate_reduction_pct.to_numpy(float)\n diffs = np.diff(vals)\n if len(diffs) and np.all(diffs <= 0):\n trend = \"nonincreasing\"\n elif len(diffs) and np.all(diffs >= 0):\n trend = \"nondecreasing\"\n else:\n trend = \"nonmonotonic\"\n local_minima = []\n for i in range(1, len(q)-1):\n if vals[i] < vals[i-1] and vals[i] < vals[i+1]:\n local_minima.append({\n \"interval\": str(q.iloc[i].interval_label),\n \"mean\": float(vals[i]),\n \"depth_vs_adjacent_mean\": float(0.5*(vals[i-1]+vals[i+1]) - vals[i]),\n })\n imin = int(np.argmin(vals)) if len(vals) else -1\n imax = int(np.argmax(vals)) if len(vals) else -1\n shape_rows.append({\n \"dataset\": ds, \"compressor\": codec, \"intervals_analyzed\": int(len(q)),\n \"trend\": trend,\n \"minimum_interval\": str(q.iloc[imin].interval_label) if imin >= 0 else \"\",\n \"minimum_mean_rate_reduction_pct\": float(vals[imin]) if imin >= 0 else np.nan,\n \"maximum_interval\": str(q.iloc[imax].interval_label) if imax >= 0 else \"\",\n \"maximum_mean_rate_reduction_pct\": float(vals[imax]) if imax >= 0 else np.nan,\n \"interior_local_minima_count\": len(local_minima),\n \"interior_local_minima\": json.dumps(local_minima, sort_keys=True),\n })\n return pdf, sdf, pd.DataFrame(shape_rows)\n\n\ndef _plot_measured_psnr_intervals(summary: pd.DataFrame) -> None:\n if summary is None or not len(summary):\n return\n try:\n import matplotlib.pyplot as plt\n datasets = sorted(summary.dataset.astype(str).unique())\n codecs = list(FINAL_CODECS)\n fig, axes = plt.subplots(len(datasets), 1, figsize=(8.0, 2.35*len(datasets)), squeeze=False)\n for ax, ds in zip(axes[:,0], datasets):\n for codec in codecs:\n q = summary[(summary.dataset.astype(str)==ds) & (summary.compressor.astype(str)==codec)].sort_values(\"interval_index\")\n if not len(q):\n continue\n x = 0.5*(q.interval_lo_db.to_numpy(float)+q.interval_hi_db.to_numpy(float))\n y = q.mean_rate_reduction_pct.to_numpy(float)\n lo = y-q.bootstrap95_lo.to_numpy(float)\n hi = q.bootstrap95_hi.to_numpy(float)-y\n ax.errorbar(x, y, yerr=np.vstack([lo,hi]), marker='o', linewidth=1.2, capsize=2.5, label=codec)\n ax.axhline(0.0, linewidth=0.8)\n ax.set_title(ds)\n ax.set_ylabel(\"Bit-rate reduction (%)\")\n ax.grid(True, alpha=0.2)\n axes[-1,0].set_xlabel(\"Measured PSNR (dB), interval midpoint\")\n handles, labels = axes[0,0].get_legend_handles_labels()\n if handles:\n fig.legend(handles, labels, loc=\"upper center\", ncol=len(labels), frameon=False)\n fig.tight_layout(rect=[0,0,1,0.97])\n fig.savefig(OUT/\"measured_psnr_interval_all_datasets.png\", dpi=220, bbox_inches=\"tight\")\n plt.close(fig)\n except Exception as e:\n log(f\"PSNR interval plot skipped: {type(e).__name__}: {e}\")\n\ndef _successful_raw_df(tables:Dict[str,List[dict]]) -> pd.DataFrame:\n raw=pd.DataFrame(tables.get(\"raw_rd_points\", []))\n if not len(raw) or \"bits_per_value\" not in raw.columns:\n return pd.DataFrame()\n raw=raw.copy()\n raw[\"bits_per_value\"]=pd.to_numeric(raw[\"bits_per_value\"],errors=\"coerce\")\n if \"psnr_db\" in raw.columns: raw[\"psnr_db\"]=pd.to_numeric(raw[\"psnr_db\"],errors=\"coerce\")\n if \"bound_ratio\" in raw.columns: raw[\"bound_ratio\"]=pd.to_numeric(raw[\"bound_ratio\"],errors=\"coerce\")\n return raw[raw[\"bits_per_value\"].notna()].copy()\n\n\ndef preliminary_codec_unit(tables:Dict[str,List[dict]],dataset:str,unit:str,codec:str) -> Optional[dict]:\n raw=_successful_raw_df(tables)\n if not len(raw): return None\n q=raw[(raw.dataset.astype(str)==str(dataset)) & (raw.unit.astype(str)==str(unit)) & (raw.compressor.astype(str)==str(codec))].copy()\n if not len(q): return None\n ref=q[q.layout_kind==\"reference\"].copy(); sel=q[q.layout_kind==\"selected\"].copy()\n out={\"reference_points\":len(ref),\"selected_points\":len(sel)}\n if len(ref) and len(sel):\n a=sel[[\"nominal_psnr\",\"bits_per_value\"]].rename(columns={\"bits_per_value\":\"sel_bpv\"})\n b=ref[[\"nominal_psnr\",\"bits_per_value\"]].rename(columns={\"bits_per_value\":\"ref_bpv\"})\n m=a.merge(b,on=\"nominal_psnr\",how=\"inner\")\n if len(m):\n sv=100.0*(1.0-m.sel_bpv.to_numpy(float)/m.ref_bpv.to_numpy(float))\n out.update({\"matched_tolerance_points\":len(m),\"same_tolerance_mean_saving_pct\":float(np.mean(sv)),\n \"same_tolerance_min_saving_pct\":float(np.min(sv)),\"same_tolerance_max_saving_pct\":float(np.max(sv))})\n integ=integrated_saving(sel,ref)\n if integ is not None: out[\"integrated_rate_saving_pct\"]=float(integ)\n if \"bound_ratio\" in q.columns and q.bound_ratio.notna().any(): out[\"max_bound_ratio\"]=float(q.bound_ratio.max())\n return out\n\n\ndef print_preliminary_codec_unit(tables:Dict[str,List[dict]],dataset:str,unit:str,codec:str) -> None:\n r=preliminary_codec_unit(tables,dataset,unit,codec)\n if not r: return\n parts=[f\"points ref/sel={r.get('reference_points',0)}/{r.get('selected_points',0)}\"]\n if \"same_tolerance_mean_saving_pct\" in r:\n parts.append(f\"same-tolerance saving mean={r['same_tolerance_mean_saving_pct']:+.3f}% min={r['same_tolerance_min_saving_pct']:+.3f}% max={r['same_tolerance_max_saving_pct']:+.3f}%\")\n if \"integrated_rate_saving_pct\" in r:\n parts.append(f\"preliminary integrated RD saving={r['integrated_rate_saving_pct']:+.3f}%\")\n if \"max_bound_ratio\" in r: parts.append(f\"max bound ratio={r['max_bound_ratio']:.8f}\")\n progress_message(\"PRELIMINARY\",f\"{dataset}/{unit}/{codec}\",\"; \".join(parts))\n\n\ndef print_preliminary_unit_table(tables:Dict[str,List[dict]],dataset:str,unit:str,codecs:Sequence[str]) -> None:\n lines=[]\n for codec in codecs:\n r=preliminary_codec_unit(tables,dataset,unit,codec)\n if not r: continue\n integ=r.get(\"integrated_rate_saving_pct\",float(\"nan\")); same=r.get(\"same_tolerance_mean_saving_pct\",float(\"nan\")); br=r.get(\"max_bound_ratio\",float(\"nan\"))\n lines.append(f\"{codec:8s} ref/sel={r.get('reference_points',0):2d}/{r.get('selected_points',0):2d} | same-tol={same:+8.3f}% | integrated={integ:+8.3f}% | max-bound={br:.8f}\")\n if lines:\n progress_message(\"UNIT SUMMARY\",f\"Preliminary results for {dataset}/{unit}\",\"\\n \".join(lines))\n\n\ndef print_preliminary_dataset_table(tables:Dict[str,List[dict]],dataset:str,codecs:Sequence[str]) -> None:\n raw=_successful_raw_df(tables)\n if not len(raw): return\n lines=[]\n units=sorted(set(raw.loc[raw.dataset.astype(str)==str(dataset),\"unit\"].astype(str)))\n for codec in codecs:\n vals=[]\n for unit in units:\n r=preliminary_codec_unit(tables,dataset,unit,codec)\n if r and \"integrated_rate_saving_pct\" in r: vals.append(r[\"integrated_rate_saving_pct\"])\n if vals:\n pos=sum(float(v)>0 for v in vals)\n lines.append(f\"{codec:8s} completed-units={len(vals):2d} | positive={pos}/{len(vals)} | mean integrated saving={np.mean(vals):+8.3f}% | median={np.median(vals):+8.3f}% | worst={np.min(vals):+8.3f}% | best={np.max(vals):+8.3f}%\")\n if lines:\n progress_message(\"DATASET SUMMARY\",f\"Running held-out summary for {dataset}\",\"\\n \".join(lines))\n\n\ndef print_selected_layout_overview(bundle:DatasetBundle, selected_by_codec:Dict[str,Layout], ranked_by_codec:Dict[str,List[Layout]]) -> None:\n lines=[]\n for codec,lay in selected_by_codec.items():\n ranks=ranked_by_codec.get(codec,[])\n top=\", \".join(x.label(bundle.axis_names) for x in ranks[:3]) if ranks else \"n/a\"\n lines.append(f\"{codec:8s} selected={lay.label(bundle.axis_names)} | top validation layouts: {top}\")\n if lines:\n progress_message(\"SELECTED LAYOUTS\",f\"{bundle.dataset_id}: validation-only decisions\",\"\\n \".join(lines))\n\n# -----------------------------------------------------------------------------\n# Benchmark execution\n# -----------------------------------------------------------------------------\ndef evaluate_dataset(bundle:DatasetBundle,adapters:List[Adapter],tables:Dict[str,List[dict]])->None:\n dr=development_range(bundle); log(f\"\\n=== {bundle.dataset_name}: shape search={bundle.search.shape}, d={bundle.d}, development range={dr:.8g} ===\")\n emit_event(\"dataset\",\"START\",dataset=bundle.dataset_id,detail=f\"search shape={bundle.search.shape}; test units={len(bundle.tests)}\")\n # unit and normalization manifests\n tables[\"dataset_units\"].append({\"dataset\":bundle.dataset_id,\"unit\":\"search\",\"role\":\"search\",\"shape\":str(tuple(bundle.search.shape)),\"values\":bundle.search.size,\"source\":bundle.source_detail})\n tables[\"dataset_units\"].append({\"dataset\":bundle.dataset_id,\"unit\":\"validation\",\"role\":\"validation\",\"shape\":str(tuple(bundle.validation.shape)),\"values\":bundle.validation.size,\"source\":bundle.source_detail})\n for u in bundle.tests: tables[\"dataset_units\"].append({\"dataset\":bundle.dataset_id,\"unit\":u.unit_id,\"role\":\"test\",\"shape\":str(tuple(u.x.shape)),\"values\":u.x.size,\"source\":bundle.source_detail})\n for c,n in enumerate(bundle.channel_names): tables[\"normalization\"].append({\"dataset\":bundle.dataset_id,\"channel\":n,\"mean\":bundle.normalization_mean[c],\"std\":bundle.normalization_std[c],\"fit_role\":\"search only\",\"ignore_zero\":bundle.normalization_ignore_zero})\n checkpoint_tables(tables)\n selected_by_codec={}; ranked_by_codec={}\n planned_rd=len(bundle.tests)*len([a for a in adapters if min(a.max_order,bundle.d)>=1])*2*len(NOMINAL_PSNR_GRID)\n existing_raw=_successful_raw_df(tables)\n existing_here=0 if not len(existing_raw) else int(((existing_raw.dataset.astype(str)==bundle.dataset_id)).sum())\n progress_message(\"BENCHMARK\",f\"{bundle.dataset_id}: held-out RD plan\",f\"planned points={planned_rd}; already checkpointed successful points for this dataset={existing_here}; test units={len(bundle.tests)}; tolerances={len(NOMINAL_PSNR_GRID)}; codecs={','.join(a.name for a in adapters)}\")\n eligible_adapters=[]\n for adapter in adapters:\n if min(adapter.max_order,bundle.d)<1: continue\n try:\n preflight_ok, preflight_detail = codec_fullshape_preflight(bundle,adapter,dr)\n tables[\"codec_eligibility\"].append({\n \"dataset\":bundle.dataset_id, \"compressor\":adapter.name,\n \"compressor_version\":adapter.version(), \"execution_target\":adapter.execution_target,\n \"required\":adapter.name in REQUIRED_CODECS,\n \"eligible\":bool(preflight_ok), \"reason\":preflight_detail,\n })\n checkpoint_tables(tables)\n if not preflight_ok:\n continue\n eligible_adapters.append(adapter)\n resumed = resumed_layout(bundle, adapter, tables)\n if resumed is not None:\n sel, ranked = resumed\n emit_event(\"layout_search\",\"SKIP\",dataset=bundle.dataset_id,compressor=adapter.name,detail=\"verified compatible selected-layout checkpoint\")\n else:\n sel,ranked=validation_layouts(bundle,adapter,dr,tables[\"layout_search\"])\n except Exception as e:\n tables[\"codec_eligibility\"].append({\n \"dataset\":bundle.dataset_id, \"compressor\":adapter.name,\n \"compressor_version\":adapter.version(), \"execution_target\":adapter.execution_target,\n \"required\":adapter.name in REQUIRED_CODECS,\n \"eligible\":False, \"reason\":str(e),\n })\n emit_event(\"layout_search\",\"FAIL\",dataset=bundle.dataset_id,compressor=adapter.name,detail=str(e))\n checkpoint_tables(tables)\n if adapter.name in OPTIONAL_CODECS:\n progress_message(\"CODEC EXCLUDED\",f\"{bundle.dataset_id}/{adapter.name}: optional baseline skipped\",str(e))\n continue\n raise\n selected_by_codec[adapter.name]=sel; ranked_by_codec[adapter.name]=ranked\n prune_stale_selected_measurements(tables,bundle.dataset_id,adapter.name,sel.machine_label())\n ref=reference_layout(bundle.d,min(adapter.max_order,bundle.d))\n tables[\"selected_layouts\"].append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"compressor_version\":adapter.version(),\"execution_target\":adapter.execution_target,\"native_order\":bundle.d,\"presented_order\":min(adapter.max_order,bundle.d),\"reference_layout\":ref.machine_label(),\"reference_semantic\":ref.label(bundle.axis_names),\"selected_layout\":sel.machine_label(),\"selected_semantic\":sel.label(bundle.axis_names),\"validation_ranked_count\":len(ranked)})\n for rank_i, rank_lay in enumerate(ranked, 1):\n tables[\"layout_rankings\"].append({\"dataset\":bundle.dataset_id,\"compressor\":adapter.name,\"selection_revision\":SELECTION_REVISION,\"compressor_version\":adapter.version(),\"execution_target\":adapter.execution_target,\"rank\":rank_i,\"layout\":rank_lay.machine_label(),\"layout_semantic\":rank_lay.label(bundle.axis_names)})\n checkpoint_tables(tables)\n print_selected_layout_overview(bundle,selected_by_codec,ranked_by_codec)\n active_names=[a.name for a in eligible_adapters if a.name in selected_by_codec]\n excluded_names=[a.name for a in adapters if a.name not in active_names and min(a.max_order,bundle.d)>=1]\n adjusted_points=len(bundle.tests)*len(active_names)*2*len(NOMINAL_PSNR_GRID)\n progress_message(\n \"BENCHMARK PLAN LOCKED\", f\"{bundle.dataset_id}: eligible codecs determined\",\n f\"active={','.join(active_names) if active_names else 'none'}; \"\n f\"excluded={','.join(excluded_names) if excluded_names else 'none'}; \"\n f\"adjusted held-out RD points={adjusted_points}\"\n )\n for unit_index,unit in enumerate(bundle.tests,1):\n emit_event(\"test_unit\",\"START\",dataset=bundle.dataset_id,unit=unit.unit_id)\n progress_message(\"TEST UNIT\",f\"{bundle.dataset_id}: starting {unit.unit_id} ({unit_index}/{len(bundle.tests)})\",f\"shape={tuple(unit.x.shape)}; values={unit.x.size:,}; each codec evaluates reference and selected layouts over {len(NOMINAL_PSNR_GRID)} tolerances\")\n for adapter in eligible_adapters:\n if adapter.name not in selected_by_codec: continue\n m=min(adapter.max_order,bundle.d); ref=reference_layout(bundle.d,m); opt=selected_by_codec[adapter.name]\n progress_message(\"CODEC\",f\"{bundle.dataset_id}/{unit.unit_id}/{adapter.name}: starting held-out sweep\",f\"reference={ref.label(bundle.axis_names)}; selected={opt.label(bundle.axis_names)}; presented order={m}\")\n geom_keys={(str(r.get(\"dataset\")),str(r.get(\"unit\")),str(r.get(\"compressor\")),str(r.get(\"layout_kind\"))) for r in tables.get(\"geometry_metrics\",[])}\n if (bundle.dataset_id,unit.unit_id,adapter.name,\"reference\") not in geom_keys:\n tables[\"geometry_metrics\"].extend(geometry_rows(bundle,unit,ref,\"reference\",adapter.name))\n if (bundle.dataset_id,unit.unit_id,adapter.name,\"selected\") not in geom_keys:\n tables[\"geometry_metrics\"].extend(geometry_rows(bundle,unit,opt,\"selected\",adapter.name))\n work=ROOT/\"tmp\"/bundle.dataset_id/adapter.name.replace(\"/\",\"_\"); work.mkdir(parents=True,exist_ok=True)\n for kind,lay in ((\"reference\",ref),(\"selected\",opt)):\n for db in NOMINAL_PSNR_GRID:\n eps=eps_from_nominal(db,dr)\n base={\"dataset\":bundle.dataset_id,\"unit\":unit.unit_id,\"compressor\":adapter.name,\"compressor_version\":adapter.version(),\"execution_target\":adapter.execution_target,\"layout_kind\":kind,\"layout\":lay.machine_label(),\"layout_semantic\":lay.label(bundle.axis_names),\"native_order\":bundle.d,\"presented_order\":m,\"nominal_psnr\":db,\"absolute_tolerance\":eps,\"development_range\":dr}\n if successful_rd_key_exists(tables,bundle.dataset_id,unit.unit_id,adapter.name,kind,float(db),lay.machine_label(),adapter.version()):\n emit_event(\"rd_point\",\"SKIP\",dataset=bundle.dataset_id,compressor=adapter.name,unit=unit.unit_id,detail=f\"{kind}; nominal={db:g}; verified successful checkpoint\")\n continue\n try:\n row,rec=one_point(adapter,unit.x,lay,eps,dr,work); full={**base,**row}; tables[\"raw_rd_points\"].append(full); tables[\"per_channel_metrics\"].extend(channel_metric_rows(bundle,unit,rec,base,kind))\n checkpoint_tables(tables)\n emit_event(\"rd_point\",\"OK\",dataset=bundle.dataset_id,compressor=adapter.name,unit=unit.unit_id,detail=f\"{kind}; nominal={db:g}; bpv={row['bits_per_value']:.6g}; measured PSNR={row['psnr_db']:.3f}; bound ratio={row['bound_ratio']:.8f}\")\n raw_now=_successful_raw_df(tables)\n ds_done=int(((raw_now.dataset.astype(str)==bundle.dataset_id)).sum()) if len(raw_now) else 0\n extra=\"\"\n if kind==\"selected\":\n refq=raw_now[(raw_now.dataset.astype(str)==bundle.dataset_id)&(raw_now.unit.astype(str)==unit.unit_id)&(raw_now.compressor.astype(str)==adapter.name)&(raw_now.layout_kind.astype(str)==\"reference\")&(pd.to_numeric(raw_now.nominal_psnr,errors=\"coerce\")==float(db))]\n if len(refq):\n ref_bpv=float(refq.iloc[-1].bits_per_value); extra=f\"; same-tolerance saving={100.0*(1.0-row['bits_per_value']/ref_bpv):+.3f}%\"\n pct=100.0*ds_done/max(1,planned_rd)\n progress_message(\"RD POINT\",f\"{bundle.dataset_id}/{unit.unit_id}/{adapter.name} {kind} nominal={db:g}\",f\"dataset progress={ds_done}/{planned_rd} ({pct:.1f}%) successful points; bpv={row['bits_per_value']:.6f}; measured PSNR={row['psnr_db']:.3f} dB; maxerr={row['max_abs_error']:.6g}; bound={row['bound_ratio']:.8f}; encode={row['encode_seconds']:.3f}s; decode={row['decode_seconds']:.3f}s{extra}\")\n if kind==\"selected\":\n print_preliminary_codec_unit(tables,bundle.dataset_id,unit.unit_id,adapter.name)\n except Exception as e:\n full={**base,\"error\":str(e)[:1000]}; tables[\"raw_rd_points\"].append(full); checkpoint_tables(tables); log(f\"ERROR {bundle.dataset_id}/{unit.unit_id}/{adapter.name}/{kind}/{db}: {e}\")\n emit_event(\"rd_point\",\"FAIL\",dataset=bundle.dataset_id,compressor=adapter.name,unit=unit.unit_id,detail=f\"{kind}; nominal={db:g}; {e}\")\n print_preliminary_codec_unit(tables,bundle.dataset_id,unit.unit_id,adapter.name)\n _plot_progress(tables,bundle.dataset_id)\n print_preliminary_unit_table(tables,bundle.dataset_id,unit.unit_id,[a.name for a in eligible_adapters if a.name in selected_by_codec])\n checkpoint_tables(tables)\n print_preliminary_dataset_table(tables,bundle.dataset_id,[a.name for a in eligible_adapters if a.name in selected_by_codec])\n emit_event(\"test_unit\",\"OK\",dataset=bundle.dataset_id,unit=unit.unit_id)\n _plot_progress(tables,bundle.dataset_id)\n emit_event(\"dataset\",\"OK\",dataset=bundle.dataset_id,detail=f\"completed {len(bundle.tests)} held-out units\")\n\n\ndef postprocess(tables:Dict[str,List[dict]])->None:\n raw=pd.DataFrame(tables[\"raw_rd_points\"])\n def _infer_target(r):\n v=str(r.get(\"execution_target\",\"\")).strip().lower()\n if v == \"cpu\": return v\n return \"cpu\"\n if len(raw):\n if \"execution_target\" not in raw.columns:\n raw[\"execution_target\"] = np.nan\n raw[\"execution_target\"] = raw.apply(_infer_target,axis=1)\n raw.to_csv(OUT/\"raw_rd_points.csv\",index=False)\n valid=raw[raw.get(\"bits_per_value\",pd.Series(index=raw.index,dtype=float)).notna()].copy()\n unit_rows=[]\n for (ds,unit,codec),g in valid.groupby([\"dataset\",\"unit\",\"compressor\"]):\n a=g[g.layout_kind==\"selected\"]; b=g[g.layout_kind==\"reference\"]; s=integrated_saving(a,b)\n if s is not None: unit_rows.append({\"dataset\":ds,\"unit\":unit,\"compressor\":codec,\"integrated_log_rate_saving_pct\":s,\"selected_psnr_lo\":a.psnr_db.min(),\"selected_psnr_hi\":a.psnr_db.max(),\"reference_psnr_lo\":b.psnr_db.min(),\"reference_psnr_hi\":b.psnr_db.max()})\n udf=pd.DataFrame(unit_rows); udf.to_csv(OUT/\"per_unit_rate_savings.csv\",index=False)\n sums=[]\n if len(udf):\n for (ds,codec),g in udf.groupby([\"dataset\",\"compressor\"]):\n x=g.integrated_log_rate_saving_pct.to_numpy(float); lo,hi=bootstrap_ci(x)\n sums.append({\"dataset\":ds,\"compressor\":codec,\"test_units\":len(x),\"mean_saving_pct\":np.mean(x),\"median_saving_pct\":np.median(x),\"min_saving_pct\":np.min(x),\"max_saving_pct\":np.max(x),\"std_saving_pct\":np.std(x,ddof=1) if len(x)>1 else 0.0,\"bootstrap95_lo\":lo,\"bootstrap95_hi\":hi,\"positive_units\":int(np.sum(x>0))})\n sdf=pd.DataFrame(sums); sdf.to_csv(OUT/\"rate_saving_summary.csv\",index=False)\n # Fresh run: all domains are produced under one protocol, so this is already the combined table.\n if len(sdf):\n sdf.assign(evidence_source=\"fresh_uniform_fast_protocol\").to_csv(OUT/\"combined_rate_saving_summary.csv\",index=False)\n if len(sdf):\n lines=[]\n for _,r in sdf.sort_values([\"dataset\",\"compressor\"]).iterrows():\n lines.append(f\"{str(r['dataset']):18s} {str(r['compressor']):8s} n={int(r['test_units'])} | mean={float(r['mean_saving_pct']):+.3f}% | median={float(r['median_saving_pct']):+.3f}% | worst={float(r['min_saving_pct']):+.3f}% | 95% CI=[{float(r['bootstrap95_lo']):+.3f},{float(r['bootstrap95_hi']):+.3f}]% | positive={int(r['positive_units'])}/{int(r['test_units'])}\")\n progress_message(\"FINAL SCIENTIFIC SUMMARY\",\"Integrated matched-PSNR rate savings\",\"\\n \".join(lines))\n\n interval_unit, interval_summary, interval_shape = measured_psnr_interval_tables(valid)\n interval_unit.to_csv(OUT/\"measured_psnr_interval_per_unit.csv\", index=False)\n interval_summary.to_csv(OUT/\"measured_psnr_interval_summary.csv\", index=False)\n interval_shape.to_csv(OUT/\"measured_psnr_interval_shape.csv\", index=False)\n _plot_measured_psnr_intervals(interval_summary)\n if len(interval_summary):\n lines=[]\n for (ds,codec), g in interval_summary.groupby([\"dataset\",\"compressor\"], sort=True):\n q=g.sort_values(\"interval_index\")\n parts=[f\"{r.interval_label}: {float(r.mean_rate_reduction_pct):+.3f}% (n={int(r.test_units)})\" for _,r in q.iterrows()]\n lines.append(f\"{ds}/{codec}: \" + \"; \".join(parts))\n progress_message(\"MEASURED-PSNR INTERVAL SUMMARY\",\"Mean bit-rate reduction over independent held-out units\",\"\\n \".join(lines))\n audits=[]\n if len(valid):\n for (ds,codec,kind),g in valid.groupby([\"dataset\",\"compressor\",\"layout_kind\"]): audits.append({\"dataset\":ds,\"compressor\":codec,\"layout_kind\":kind,\"points\":len(g),\"max_bound_ratio\":g.bound_ratio.max(),\"mean_bound_ratio\":g.bound_ratio.mean(),\"violations_gt_1_plus_1e6\":int((g.bound_ratio>1.000001).sum())})\n pd.DataFrame(audits).to_csv(OUT/\"error_bound_audit.csv\",index=False)\n rt=[]\n if len(valid):\n for (ds,codec,kind),g in valid.groupby([\"dataset\",\"compressor\",\"layout_kind\"]):\n target = str(g.execution_target.iloc[0]) if \"execution_target\" in g.columns and len(g) else \"unknown\"\n rt.append({\"dataset\":ds,\"compressor\":codec,\"execution_target\":target,\"layout_kind\":kind,\"points\":len(g),\"median_encode_seconds\":g.encode_seconds.median(),\"median_decode_seconds\":g.decode_seconds.median(),\"median_encode_MBps\":g.encode_MBps.median(),\"median_decode_MBps\":g.decode_MBps.median()})\n pd.DataFrame(rt).to_csv(OUT/\"runtime_summary.csv\",index=False)\n\n # Cross-compressor comparisons use validation-selected single layouts only.\n cross=[]; envrows=[]\n if len(valid):\n for (ds,unit),gu in valid[valid.layout_kind==\"selected\"].groupby([\"dataset\",\"unit\"]):\n curves={c:g.sort_values(\"psnr_db\") for c,g in gu.groupby(\"compressor\") if len(g)>=2}\n if \"SZ3\" in curves:\n sz=curves[\"SZ3\"]\n for c,g in curves.items():\n if c==\"SZ3\": continue\n sv=integrated_saving(g,sz) # positive => contender uses fewer bits than selected SZ3\n if sv is not None: cross.append({\"dataset\":ds,\"unit\":unit,\"compressor\":c,\"integrated_rate_saving_vs_selected_SZ3_pct\":sv})\n ext={c:g for c,g in curves.items() if c!=\"SZ3\"}\n if ext:\n lo=float(sz.psnr_db.min()); hi=float(sz.psnr_db.max()); targets=np.linspace(lo,hi,101); lr=[]; used=0\n for q in targets:\n if qsz.psnr_db.max(): continue\n rsz=float(np.interp(q,sz.psnr_db,sz.bits_per_value)); opts=[]\n for c,g in ext.items():\n if q>=g.psnr_db.min() and q<=g.psnr_db.max(): opts.append(float(np.interp(q,g.psnr_db,g.bits_per_value)))\n if opts:\n rext=min(opts); lr.append(math.log(max(rsz,1e-30))-math.log(max(rext,1e-30))); used+=1\n if lr:\n # positive means selected SZ3 saves rate relative to the evaluated external envelope\n envrows.append({\"dataset\":ds,\"unit\":unit,\"matched_grid_points\":used,\"selected_SZ3_saving_vs_external_envelope_pct\":100.0*(1.0-math.exp(float(np.mean(lr))))})\n pd.DataFrame(cross).to_csv(OUT/\"cross_codec_pairwise.csv\",index=False)\n pd.DataFrame(envrows).to_csv(OUT/\"selected_sz3_vs_external_envelope.csv\",index=False)\n\n # Other direct tables\n for name in (\"dataset_units\",\"normalization\",\"layout_search\",\"layout_rankings\",\"selected_layouts\",\"per_channel_metrics\",\"geometry_metrics\",\"codec_eligibility\"):\n pd.DataFrame(tables[name]).to_csv(OUT/(name+\".csv\"),index=False)\n _plot_progress(tables)\n _plot_final_summaries()\n emit_event(\"postprocess\",\"DONE\",detail=\"summary tables and diagnostic plots written\")\n\n\ndef environment_record() -> dict:\n d={\"timestamp_utc\":time.strftime(\"%Y-%m-%dT%H:%M:%SZ\",time.gmtime()),\"python\":sys.version,\"platform\":platform.platform(),\"machine\":platform.machine(),\"processor\":platform.processor(),\"seed\":SEED,\"cpu_count\":os.cpu_count(),\"packages\":{n:package_version(n) for n in (\"numpy\",\"pandas\",\"pysz\",\"hdf5plugin\",\"h5py\",\"nibabel\",\"netCDF4\")}}\n try:\n import psutil; d[\"ram_gb\"]=psutil.virtual_memory().total/2**30\n except Exception: pass\n d[\"execution_scope\"] = \"five-domain two-codec CPU benchmark: SZ3,HPEZ\"\n return d\n\n\ndef write_protocol() -> None:\n protocol={\n \"primary_question\":\"rate change caused by validation-selected dimension permutation and contiguous folding around an unchanged compressor\",\n \"compressors\":{\"required\":[\"SZ3\",\"HPEZ\"],\"optional\":[]},\n \"claim_scope\":\"representation selection across two prediction/interpolation-based error-bounded compressors; no claim of universality across all codec architectures\",\n \"execution_scope\":\"uniform CPU-only protocol across five required domains; no dataset-specific non-reference layout is injected\",\n \"datasets\":DATASET_URLS,\n \"fresh_protocol\":True,\n \"headline_selection\":\"one validation-selected layout per dataset and compressor; no held-out selection\",\n \"reference_layout\":\"native semantic axis order; fold the first d-m+1 axes and keep remaining axes explicit\",\n \"layout_space\":\"all semantic-axis permutations times all contiguous partitions into the compressor-supported order\",\n \"proxy_policy\":\"exhaustive locality proxy shared once per dataset/interface-order; codec-specific coarse and full-validation probes remain separate\",\n \"proxy_global_keep\":PROXY_GLOBAL_KEEP,\"proxy_per_partition_keep\":PROXY_PER_PARTITION_KEEP,\n \"coarse_keep\":COARSE_KEEP,\"final_keep\":FINAL_KEEP,\n \"selection_revision\":SELECTION_REVISION,\"validation_winner_grid\":VALID_WINNER_GRID,\n \"quality_grid_nominal_psnr\":NOMINAL_PSNR_GRID.tolist(),\n \"search_grid_nominal_psnr\":SEARCH_NOMINAL_GRID.tolist(),\n \"validation_grid_nominal_psnr\":VALID_NOMINAL_GRID.tolist(),\n \"quality_mapping\":\"eps = development_range * 10^(-nominal_psnr/20)\",\n \"reported_quality\":\"independently measured PSNR/RMSE/NRMSE and maximum absolute error\",\n \"error_contract\":f\"retain only max_abs_error / requested_abs_tolerance <= {BOUND_ACCEPT_RATIO}\",\n \"rate\":\"8*(codec_payload + common_wrapper_bytes)/number_of_values\",\n \"common_wrapper_bytes\":COMMON_WRAPPER_BYTES,\n \"statistical_unit\":\"non-overlapping held-out temporal block, held-out subject, or held-out plot\",\n \"summary_metric\":\"integrated mean log-rate ratio over common measured-PSNR interval; converted to percent bitrate saving\",\n \"candidate_policy\":\"the reference layout is retained as the baseline; every non-reference candidate is produced uniformly by exhaustive enumeration, global plus per-fold-partition locality screening, codec-specific coarse screening, and full validation for every dataset and compressor\",\n \"proxy_shortlist_policy\":\"union of the globally lowest locality scores and the lowest locality scores within every contiguous fold partition; no semantic layout identity is specified\",\n \"dataset_specific_forced_nonreference_candidates\":[],\n \"build_parallel_jobs\":BUILD_JOBS,\n \"execution\":\"CPU-only; no GPU build\",\n \"memory_policy\":\"one dataset resident at a time; bounded error-metric temporaries; deterministic probe crops\",\n \"resume_policy\":\"checkpoints are written after selection and every RD point; --fresh/HOC_FRESH=1 deletes all scientific state but preserves the separately verified codec-binary cache\",\n \"measured_psnr_interval_analysis\":{\"start_db\":PSNR_INTERVAL_START_DB,\"width_db\":PSNR_INTERVAL_WIDTH_DB,\"count\":PSNR_INTERVAL_COUNT,\"grid_points\":PSNR_INTERVAL_GRID_POINTS,\"minimum_coverage_db\":PSNR_INTERVAL_MIN_COVERAGE_DB,\"role\":\"held-out descriptive analysis only; never used for layout selection\"},\n \"checkpoint_schema\":CHECKPOINT_SCHEMA\n }\n (OUT/\"protocol.json\").write_text(json.dumps(protocol,indent=2),encoding=\"utf-8\")\n\ndef write_selection_integrity() -> None:\n audit = {\n \"selection_revision\": SELECTION_REVISION,\n \"datasets\": [\"ams_gefs\", \"mrms\", \"instant_odc\", \"hyperleaf\", \"ixi\"],\n \"compressors\": list(FINAL_CODECS),\n \"reference_layout_retained_as_baseline\": True,\n \"dataset_specific_forced_nonreference_candidates\": [],\n \"nonreference_candidate_source\": \"uniform exhaustive permutation/folding enumeration -> locality screen -> codec coarse screen -> full validation selection\",\n \"held_out_measurements_used_for_selection\": False,\n \"selection_rule_identical_across_dataset_codec_pairs\": True,\n \"note\": \"Protocol constants are fixed before execution; no semantic layout identity is specified anywhere in selection. Layout identities and rankings are data-dependent.\"\n }\n (OUT / \"selection_integrity.json\").write_text(json.dumps(audit, indent=2), encoding=\"utf-8\")\n\n\ndef evidence_zip() -> Path:\n zpath=ROOT/\"high_order_compression_evidence.zip\"\n keep=[p for p in OUT.iterdir() if p.is_file() and p.suffix.lower() in {\".csv\",\".json\",\".txt\",\".png\"}]\n with zipfile.ZipFile(zpath,\"w\",compression=zipfile.ZIP_DEFLATED,compresslevel=9) as zf:\n for p in sorted(keep): zf.write(p,arcname=p.name)\n return zpath\n\n\n\ndef verify_all_codecs_before_datasets(adapters: Sequence[Adapter]) -> pd.DataFrame:\n \"\"\"Strict executable+numerical gate before any real dataset is loaded.\n\n This is intentionally stronger than checking imports or process exit codes.\n Every configured codec must round-trip deterministic, asymmetric float32\n tensors at three absolute tolerances and satisfy the common pointwise bound.\n The asymmetric construction makes axis/order mistakes fail conspicuously.\n \"\"\"\n expected = tuple(FINAL_CODECS)\n by_name = {a.name: a for a in adapters}\n rows = []\n\n progress_message(\n \"CODEC VALIDATION GATE\",\n \"Verifying all CPU codecs numerically before dataset work\",\n \"Required for this notebook run: \" + \", \".join(expected) +\n \". Each codec must pass 3 absolute-error round trips; otherwise the benchmark stops before loading data.\"\n )\n\n log(\"[CPU CHECK] CPU-only codec gate; no GPU dependency\")\n\n missing = [name for name in expected if name not in by_name]\n if missing:\n for name in missing:\n rows.append({\n \"compressor\": name, \"trial\": \"setup\", \"epsilon\": np.nan,\n \"payload_bytes\": np.nan, \"bits_per_value\": np.nan,\n \"max_abs_error\": np.nan, \"bound_ratio\": np.inf,\n \"shape_ok\": False, \"finite_ok\": False, \"status\": \"FAIL\",\n \"detail\": \"codec unavailable after setup/build\"\n })\n pd.DataFrame(rows).to_csv(OUT/\"codec_validation_gate.csv\", index=False)\n status = {}\n try:\n sdf = pd.read_csv(OUT/\"software_status.csv\")\n for _, rr in sdf.iterrows():\n status[str(rr.get(\"compressor\"))] = str(rr.get(\"detail\", \"\"))\n except Exception:\n pass\n details = \"; \".join(f\"{name}: {status.get(name, 'setup/build failed')}\" for name in missing)\n raise RuntimeError(\n \"STRICT CPU CODEC GATE FAILED: unavailable codecs: \" + \", \".join(missing) +\n f\". Details: {details}. No GPU is required for this benchmark.\"\n )\n\n workroot = ROOT / \"codec_validation_gate\"\n shutil.rmtree(workroot, ignore_errors=True)\n workroot.mkdir(parents=True, exist_ok=True)\n\n # Do not use symmetric shapes or a simple linear ramp: an ordering bug can\n # otherwise survive a weak round-trip test.\n for name in expected:\n adapter = by_name[name]\n shape = (4, 7, 11, 13) if adapter.max_order >= 4 else (17, 19, 23)\n idx = np.indices(shape, dtype=np.float32)\n a = np.zeros(shape, dtype=np.float32)\n for j, g in enumerate(idx):\n a += (j + 1.37) * g\n a += np.sin(g * (0.071 + 0.019*j)).astype(np.float32) * (0.31 + 0.11*j)\n # Add coordinate-coupled structure and deterministic sparse impulses.\n flat = a.reshape(-1)\n flat[::37] += 3.25\n flat[11::101] -= 1.75\n a = np.ascontiguousarray(a, dtype=np.float32)\n dr = float(np.max(a) - np.min(a))\n eps_list = [dr * 1e-2, dr * 2e-3, dr * 5e-4]\n\n log(f\"[CODEC GATE] {name}: version={adapter.version()} target={adapter.execution_target} shape={shape} range={dr:.9g}\")\n codec_ok = True\n for trial, eps in enumerate(eps_list, 1):\n sub = workroot / name / f\"trial_{trial}\"\n sub.mkdir(parents=True, exist_ok=True)\n try:\n r = adapter.compress(a, float(eps), sub)\n rec = np.asarray(r.reconstruction, dtype=np.float32)\n shape_ok = tuple(rec.shape) == tuple(a.shape)\n finite_ok = bool(np.isfinite(rec).all())\n if shape_ok and finite_ok:\n ma = float(np.max(np.abs(a.astype(np.float64) - rec.astype(np.float64))))\n ratio = ma / float(eps) if eps > 0 else np.inf\n else:\n ma, ratio = np.inf, np.inf\n payload = int(r.payload_bytes)\n bpv = 8.0 * payload / a.size\n passed = bool(shape_ok and finite_ok and payload > 0 and np.isfinite(ratio) and ratio <= BOUND_ACCEPT_RATIO)\n rows.append({\n \"compressor\": name, \"trial\": trial, \"epsilon\": float(eps),\n \"payload_bytes\": payload, \"bits_per_value\": bpv,\n \"max_abs_error\": ma, \"bound_ratio\": ratio,\n \"shape_ok\": shape_ok, \"finite_ok\": finite_ok,\n \"status\": \"PASS\" if passed else \"FAIL\",\n \"detail\": r.detail,\n })\n log(\n f\" trial {trial}/3 eps={eps:.8g} -> bpv={bpv:.6f}; \"\n f\"maxerr={ma:.8g}; bound={ratio:.8f}; shape={'OK' if shape_ok else 'BAD'}; \"\n f\"finite={'OK' if finite_ok else 'BAD'}; {'PASS' if passed else 'FAIL'}\"\n )\n if not passed:\n codec_ok = False\n break\n except Exception as e:\n codec_ok = False\n rows.append({\n \"compressor\": name, \"trial\": trial, \"epsilon\": float(eps),\n \"payload_bytes\": np.nan, \"bits_per_value\": np.nan,\n \"max_abs_error\": np.nan, \"bound_ratio\": np.inf,\n \"shape_ok\": False, \"finite_ok\": False, \"status\": \"FAIL\",\n \"detail\": f\"{type(e).__name__}: {e}\",\n })\n log(f\" trial {trial}/3 -> FAIL: {type(e).__name__}: {e}\")\n break\n\n log(f\"[CODEC GATE RESULT] {name}: {'VALID' if codec_ok else 'INVALID'}\")\n\n df = pd.DataFrame(rows)\n df.to_csv(OUT/\"codec_validation_gate.csv\", index=False)\n good = set(df.loc[df[\"status\"].eq(\"PASS\"), \"compressor\"].astype(str))\n failed = []\n for name in expected:\n sub = df[df[\"compressor\"].eq(name)]\n if len(sub) != 3 or not sub[\"status\"].eq(\"PASS\").all():\n failed.append(name)\n\n if failed and STRICT_ALL_CODEC_GATE:\n raise RuntimeError(\n \"STRICT CODEC GATE FAILED for: \" + \", \".join(failed) +\n \". No scientific dataset has been benchmarked. Fix/enable the codec and rerun the same notebook cell. \"\n f\"Audit: {OUT/'codec_validation_gate.csv'}\"\n )\n\n progress_message(\n \"CODEC VALIDATION GATE\",\n \"BOTH CPU CODECS VALID\",\n \"SZ3 and HPEZ each passed 3/3 deterministic absolute-error round trips. \"\n f\"Audit written to {OUT/'codec_validation_gate.csv'}. Real-dataset full-shape preflight will still run before each layout search.\"\n )\n return df\n\ndef main() -> None:\n if FRESH_RUN and ROOT.exists():\n # Preserve nothing scientific: this is the explicit from-zero switch requested for final reruns.\n shutil.rmtree(ROOT,ignore_errors=True)\n for p in (ROOT,OUT,CACHE): p.mkdir(parents=True,exist_ok=True)\n TOOLS.mkdir(parents=True,exist_ok=True)\n progress_message(\"0/4 START AUDIT\",\"Preparing fresh uniform benchmark\",\n f\"fresh={FRESH_RUN}; resume={RESUME_ENABLED}; work={ROOT}; input={INPUT_ROOT}; required_for_run=SZ3,HPEZ\")\n tables,resume_state,checkpoint_compatible = load_resume_state()\n prune_nonfinal_codec_rows(tables)\n save_resume_state(resume_state)\n emit_event(\"benchmark\",\"START\",detail=f\"resume-first benchmark; compatible_checkpoint={checkpoint_compatible}\")\n progress_message(\"INITIALIZATION\",\"Benchmark started\",\n f\"work={ROOT}; input={INPUT_ROOT}; strict_datasets={STRICT_DATASETS}; strict_codecs={STRICT_CODECS}; build_jobs={BUILD_JOBS}; codecs=SZ3,HPEZ; execution=CPU-only\")\n write_protocol(); write_selection_integrity(); (OUT/\"environment.json\").write_text(json.dumps(environment_record(),indent=2),encoding=\"utf-8\")\n # Locate optional harness-side metadata without assuming a normal .py execution\n # context. Kaggle/Jupyter may execute this source with no ``__file__``.\n harness_candidates = []\n env_harness = os.environ.get(\"HOC_HARNESS_DIR\", \"\").strip()\n if env_harness:\n harness_candidates.append(Path(env_harness))\n script_name = globals().get(\"__file__\")\n if script_name:\n harness_candidates.append(Path(script_name).resolve().parent)\n harness_candidates.extend([\n Path.cwd(),\n Path(\"/kaggle/working/high_order_compression_harness\"),\n Path(\"/kaggle/working\"),\n ])\n local_manifest = next(\n (d / \"dataset_manifest.json\" for d in harness_candidates\n if (d / \"dataset_manifest.json\").exists()),\n None,\n )\n if local_manifest is not None:\n shutil.copy2(local_manifest, OUT / \"dataset_manifest.json\")\n progress_message(\"INITIALIZATION\", \"Dataset manifest found\", f\"Copied provenance manifest from {local_manifest}. Dataset discovery itself uses /kaggle/input and does not depend on this file.\")\n else:\n progress_message(\"INITIALIZATION\", \"Dataset manifest not found (non-fatal)\",\n \"The manifest is provenance/documentation only. The benchmark discovers all five configured datasets directly under /kaggle/input. Missing datasets are reported and skipped unless HOC_STRICT_DATASETS=1.\")\n paper_index={\n \"headline_representation_table\":\"rate_saving_summary.csv\",\n \"per_independent_unit_effects\":\"per_unit_rate_savings.csv\",\n \"absolute_rate_distortion_points\":\"raw_rd_points.csv\",\n \"selected_coordinate_maps\":\"selected_layouts.csv\",\n \"layout_search_audit\":\"layout_search.csv\",\n \"layout_validation_ranking\":\"layout_rankings.csv\",\n \"resume_state\":\"resume_state.json\",\n \"resume_protocol_fingerprint\":\"resume_protocol_fingerprint.json\",\n \"pointwise_error_validation\":\"error_bound_audit.csv\",\n \"systems_performance\":\"runtime_summary.csv\",\n \"scientific_variable_metrics\":\"per_channel_metrics.csv\",\n \"coordinate_geometry_mechanism\":\"geometry_metrics.csv\",\n \"cross_compressor_context\":\"cross_codec_pairwise.csv\",\n \"selected_SZ3_against_evaluated_external_envelope\":\"selected_sz3_vs_external_envelope.csv\",\n \"software_and_builds\":\"software_status.csv\",\n \"machine_environment\":\"environment.json\",\n \"exact_protocol\":\"protocol.json\",\n \"progress_events\":\"progress_events.csv\",\n \"live_completion_plot\":\"progress_completion.png\",\n \"live_rd_plot\":\"progress_rd_current.png\",\n \"live_error_bound_plot\":\"progress_error_bound.png\",\n \"combined_crossdomain_summary\":\"combined_rate_saving_summary.csv\",\n \"measured_psnr_interval_per_unit\":\"measured_psnr_interval_per_unit.csv\",\n \"measured_psnr_interval_summary\":\"measured_psnr_interval_summary.csv\",\n \"measured_psnr_interval_shape\":\"measured_psnr_interval_shape.csv\",\n \"measured_psnr_interval_plot\":\"measured_psnr_interval_all_datasets.png\",\n \"selection_integrity\":\"selection_integrity.json\",\n \"final_rate_saving_plot\":\"final_rate_savings.png\",\n \"final_throughput_plot\":\"final_encode_throughput.png\"\n }\n (OUT/\"paper_index.json\").write_text(json.dumps(paper_index,indent=2),encoding=\"utf-8\")\n adapters,_=setup_adapters()\n verify_all_codecs_before_datasets(adapters)\n verify_all_five_dataset_sources()\n progress_message(\"PRECISION PROTOCOL\",\"Expanded validation and independent-unit precision without test leakage\",\n f\"held-out RD points={NOMINAL_PSNR_GRID.tolist()}; coarse={SEARCH_NOMINAL_GRID.tolist()}; validation={VALID_NOMINAL_GRID.tolist()}; proxy_global_keep={PROXY_GLOBAL_KEEP}; proxy_per_partition_keep={PROXY_PER_PARTITION_KEEP}; coarse_keep={COARSE_KEEP}; no dataset-specific candidate; proxy shared by interface order\")\n progress_message(\"DATASETS\", \"Configured Kaggle inputs\", \"\\n \".join(f\"{k}: {v}\" for k,v in DATASET_URLS.items()))\n loaders=[(\"ams_gefs\",load_ams),(\"mrms\",load_mrms),(\"instant_odc\",load_instant_odc),(\"hyperleaf\",load_hyperleaf),(\"ixi\",load_ixi)]\n dataset_errors=[]\n # Load, evaluate, checkpoint, and release each dataset under the same protocol.\n for name,fn in loaders:\n if dataset_marked_complete(resume_state, name, adapters):\n progress_message(\"RESUME\", f\"Dataset {name} already complete — load/evaluation skipped\", \"Completion marker matches the current protocol fingerprint.\")\n emit_event(\"dataset\",\"SKIP\",dataset=name,detail=\"compatible completed-dataset checkpoint\")\n continue\n b=None\n try:\n progress_message(\"DATASET\", f\"Loading {name}\", \"Only this dataset will be resident in memory; it will be released before the next dataset is loaded.\")\n emit_event(\"dataset_load\",\"START\",dataset=name)\n b=fn(); emit_event(\"dataset_load\",\"OK\",dataset=name,detail=f\"search={b.search.shape}; validation={b.validation.shape}; tests={len(b.tests)}\")\n progress_message(\"DATASET\", f\"Loaded {name}\", f\"search={b.search.shape}; validation={b.validation.shape}; held-out units={len(b.tests)}; axes={b.axis_names}\")\n evaluate_dataset(b,adapters,tables)\n checkpoint_tables(tables)\n mark_dataset_complete(resume_state, name, adapters)\n progress_message(\"DATASET\", f\"Completed {name}\", \"Checkpoint committed. Future compatible runs will skip loading and evaluating this dataset.\")\n except Exception as e:\n dataset_errors.append({\"dataset\":name,\"error\":str(e)}); log(f\"DATASET ERROR {name}: {e}\"); emit_event(\"dataset_load\",\"FAIL\",dataset=name,detail=str(e))\n if STRICT_DATASETS:\n pd.DataFrame(dataset_errors).to_csv(OUT/\"dataset_errors.csv\",index=False)\n raise\n finally:\n if b is not None:\n del b\n gc.collect()\n pd.DataFrame(dataset_errors).to_csv(OUT/\"dataset_errors.csv\",index=False)\n checkpoint_tables(tables)\n postprocess(tables)\n # Refresh environment after builds to capture versions/status.\n (OUT/\"environment.json\").write_text(json.dumps(environment_record(),indent=2),encoding=\"utf-8\")\n emit_event(\"benchmark\",\"DONE\",detail=\"uniform five-domain benchmark completed\")\n z=evidence_zip(); log(f\"FINAL_EVIDENCE_ZIP={z}\"); print(z)\n\nif __name__==\"__main__\":\n try:\n main()\n except Exception as e:\n try:\n emit_event(\"benchmark\",\"FAIL\",detail=f\"{type(e).__name__}: {e}\")\n log(f\"FATAL ERROR: {type(e).__name__}: {e}\")\n z=evidence_zip()\n print(f\"PARTIAL_EVIDENCE_ZIP={z}\", flush=True)\n finally:\n raise\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-09-19T00:40:55.703843Z","iopub.execute_input":"2026-09-19T00:40:55.704167Z","iopub.status.idle":"2026-09-19T02:24:02.569956Z","shell.execute_reply.started":"2026-09-19T00:40:55.704136Z","shell.execute_reply":"2026-09-19T02:24:02.569003Z"}},"outputs":[],"execution_count":null}]}