diff --git a/.github/quality/workflow-inventory-budget-v1.json b/.github/quality/workflow-inventory-budget-v1.json index 01589de3d..94de84fe3 100644 --- a/.github/quality/workflow-inventory-budget-v1.json +++ b/.github/quality/workflow-inventory-budget-v1.json @@ -2,7 +2,7 @@ "schema": "bayesian-phystwin.workflow-inventory-budget", "schema_version": 1, "baseline_revision": "45e1f4454d50fd1970af13578f0383872814125e", - "maximum_checked_in_workflows": 82, + "maximum_checked_in_workflows": 83, "temporary_looking_workflow_allowlist": [], "retirement_target_maximum_checked_in_workflows": 81, "retirement_target_maximum_temporary_looking_workflows": 0 diff --git a/.github/workflows/tracking-cloth-evaluation.yml b/.github/workflows/tracking-cloth-evaluation.yml new file mode 100644 index 000000000..48bc3c871 --- /dev/null +++ b/.github/workflows/tracking-cloth-evaluation.yml @@ -0,0 +1,299 @@ +# workflow-lifecycle: permanent +# workflow-owner: IPS-Stuttgart maintainers +name: Tracking Cloth Deformation evaluation + +on: + pull_request: + paths: + - .github/workflows/tracking-cloth-evaluation.yml + - .github/requests/tracking-cloth-deformation-v1-evaluate.json + - experiments/tracking_cloth_deformation_v1/** + - tests/test_tracking_cloth_deformation_v1.py + push: + branches: [main] + paths: + - .github/requests/tracking-cloth-deformation-v1-evaluate.json + workflow_dispatch: + inputs: + mode: + description: "inventory / source-only qualification / sealed shake-to-twist pilot" + type: choice + options: [inventory, source_only, evaluate] + default: source_only + required: true + workers: + description: "Source rollout CPU workers (GPU is not needed)" + type: choice + options: ["1", "2", "4", "8"] + default: "4" + required: true + +permissions: + contents: read + +concurrency: + group: tracking-cloth-evaluation-${{ github.event_name }}-${{ github.ref }} + cancel-in-progress: false + +env: + PYTHONUNBUFFERED: "1" + PYTHONDONTWRITEBYTECODE: "1" + OPENBLAS_NUM_THREADS: "1" + OMP_NUM_THREADS: "1" + MKL_NUM_THREADS: "1" + +jobs: + contracts: + name: Tracking-cloth synthetic contracts + runs-on: ubuntu-latest + timeout-minutes: 15 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0 + with: + python-version: "3.12" + cache: pip + cache-dependency-path: experiments/tracking_cloth_deformation_v1/requirements.txt + - name: Install isolated test dependencies + run: python -m pip install -r experiments/tracking_cloth_deformation_v1/requirements.txt pytest ruff + - name: Test data boundaries, predictors, and sealed scoring + run: python -m pytest -q tests/test_tracking_cloth_deformation_v1.py + - name: Validate a canonical real-data request when present + shell: bash + run: | + python - <<'PY' + import json + from pathlib import Path + + path = Path(".github/requests/tracking-cloth-deformation-v1-evaluate.json") + if not path.exists(): + raise SystemExit(0) + request = json.loads(path.read_text()) + required = { + "schema": "tracking-cloth-deformation-evaluation-request-v1", + "mode": "evaluate", + "dataset_root": "/home/github-runner/.cache/datasets/tracking-cloth-deformation-v1-zenodo-14644526", + "runner_label": "gpuserver4090", + "workers": 4, + "authorize_target_scoring": True, + "paper_claim_authorized": False, + } + for key, value in required.items(): + if request.get(key) != value: + raise SystemExit(f"Invalid request field {key!r}: {request.get(key)!r}") + allowed = set(required) | {"request_id", "expected_source_revision", "evidence_class"} + unknown = set(request) - allowed + if unknown: + raise SystemExit(f"Unknown request fields: {sorted(unknown)}") + if not request.get("request_id"): + raise SystemExit("request_id must be nonempty") + revision = request.get("expected_source_revision", "") + if len(revision) != 40 or any(c not in "0123456789abcdef" for c in revision): + raise SystemExit("expected_source_revision must be a lowercase commit SHA") + if request.get("evidence_class") != "public-real-data-pilot; no fresh-confirmation claim": + raise SystemExit("Unexpected evidence_class") + PY + - name: Lint and formatting diagnostics + if: always() + run: | + python -m ruff format --diff experiments/tracking_cloth_deformation_v1 tests/test_tracking_cloth_deformation_v1.py > formatting.patch || true + cat formatting.patch + python -m ruff check experiments/tracking_cloth_deformation_v1 tests/test_tracking_cloth_deformation_v1.py + python -m ruff format --check experiments/tracking_cloth_deformation_v1 tests/test_tracking_cloth_deformation_v1.py + - name: Retain a formatting patch on validation failure + if: failure() + uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.0 + with: + name: tracking-cloth-formatting-${{ github.run_id }} + path: formatting.patch + retention-days: 7 + if-no-files-found: ignore + + authorize-real-data: + name: Authorize target-closed real-data request + if: github.event_name != 'pull_request' + runs-on: ubuntu-latest + timeout-minutes: 5 + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + fetch-depth: 2 + - name: Require main for manual dispatch + if: github.event_name == 'workflow_dispatch' + shell: bash + run: test "$GITHUB_REF" = refs/heads/main + - name: Require one newly added canonical request on push + if: github.event_name == 'push' + shell: bash + env: + BEFORE_SHA: ${{ github.event.before }} + PUSH_FORCED: ${{ github.event.forced }} + run: | + set -euo pipefail + request=.github/requests/tracking-cloth-deformation-v1-evaluate.json + test "$GITHUB_REF" = refs/heads/main + test "$PUSH_FORCED" = false + test "$BEFORE_SHA" != 0000000000000000000000000000000000000000 + mapfile -t changed < <(git diff --name-status "$BEFORE_SHA" "$GITHUB_SHA") + printf 'Changed files:\n%s\n' "${changed[*]}" + test "${#changed[@]}" -eq 1 + test "${changed[0]}" = $'A\t'"$request" + python - <<'PY' + import json + import os + import subprocess + from pathlib import Path + + path = Path(".github/requests/tracking-cloth-deformation-v1-evaluate.json") + request = json.loads(path.read_text()) + parent = subprocess.check_output( + ["git", "rev-parse", f"{os.environ['GITHUB_SHA']}^"], text=True + ).strip() + if request["expected_source_revision"] != parent: + raise SystemExit( + "Request expected_source_revision does not equal the trigger commit parent" + ) + PY + + evaluation: + name: Read-only cloth pilot / gpuserver4090 + needs: [contracts, authorize-real-data] + if: >- + always() && + (github.event_name == 'workflow_dispatch' || github.event_name == 'push') && + github.ref == 'refs/heads/main' && + github.repository == 'IPS-Stuttgart/BayesianPhysTwin' && + needs.contracts.result == 'success' && + needs.authorize-real-data.result == 'success' + runs-on: [self-hosted, Linux, X64, gpuserver4090] + timeout-minutes: 90 + env: + DATASET_ROOT: /home/github-runner/.cache/datasets/tracking-cloth-deformation-v1-zenodo-14644526 + EVALUATION_MODE: ${{ github.event_name == 'push' && 'evaluate' || inputs.mode }} + WORKERS: ${{ github.event_name == 'push' && '4' || inputs.workers }} + steps: + - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 + with: + persist-credentials: false + - uses: actions/setup-python@5fda3b95a4ea91299a34e894583c3862153e4b97 # v7.0.0 + with: + python-version: "3.12" + - name: Prepare private scratch and isolated environment + shell: bash + run: | + set -euo pipefail + test -d "$DATASET_ROOT" + test -r "$DATASET_ROOT" + echo "Runner: $RUNNER_NAME; required label: gpuserver4090" + echo "Included dataset license: CC BY-NC-SA 4.0; metadata conflict retained." + venv="$RUNNER_TEMP/tracking-cloth-venv-$GITHUB_RUN_ID-$GITHUB_RUN_ATTEMPT" + output="$RUNNER_TEMP/tracking-cloth-$GITHUB_RUN_ID-$GITHUB_RUN_ATTEMPT" + test ! -e "$venv" + test ! -e "$output" + python -m venv "$venv" + "$venv/bin/python" -m pip install -r experiments/tracking_cloth_deformation_v1/requirements.txt + echo "CLOTH_PY=$venv/bin/python" >> "$GITHUB_ENV" + echo "CLOTH_VENV=$venv" >> "$GITHUB_ENV" + echo "CLOTH_OUT=$output" >> "$GITHUB_ENV" + - name: Audit, fit sources, and optionally seal target predictions + shell: bash + run: | + set -euo pipefail + case "$EVALUATION_MODE" in + inventory) stage=inventory ;; + source_only) stage=source ;; + evaluate) stage=predict ;; + *) echo "Unsupported mode" >&2; exit 2 ;; + esac + "$CLOTH_PY" -m experiments.tracking_cloth_deformation_v1.run \ + --dataset-root "$DATASET_ROOT" --output "$CLOTH_OUT" \ + --stage "$stage" --workers "$WORKERS" + "$CLOTH_PY" - <<'PY' + import hashlib + import json + import os + from pathlib import Path + + event = os.environ["GITHUB_EVENT_NAME"] + request_path = Path( + ".github/requests/tracking-cloth-deformation-v1-evaluate.json" + ) + if event == "push": + raw = request_path.read_bytes() + request = json.loads(raw) + request_sha256 = hashlib.sha256(raw).hexdigest() + else: + request = { + "schema": "tracking-cloth-deformation-manual-dispatch-v1", + "request_id": f"github-run-{os.environ['GITHUB_RUN_ID']}-attempt-{os.environ['GITHUB_RUN_ATTEMPT']}", + "mode": os.environ["EVALUATION_MODE"], + "dataset_root": os.environ["DATASET_ROOT"], + "runner_label": "gpuserver4090", + "workers": int(os.environ["WORKERS"]), + "authorize_target_scoring": os.environ["EVALUATION_MODE"] == "evaluate", + "paper_claim_authorized": False, + } + raw = json.dumps(request, sort_keys=True, separators=(",", ":")).encode() + request_sha256 = hashlib.sha256(raw).hexdigest() + record = { + **request, + "request_sha256": request_sha256, + "github_sha": os.environ["GITHUB_SHA"], + "github_run_id": os.environ["GITHUB_RUN_ID"], + "github_run_attempt": os.environ["GITHUB_RUN_ATTEMPT"], + } + output = Path(os.environ["CLOTH_OUT"]) / "evaluation_request.json" + output.write_text(json.dumps(record, indent=2, sort_keys=True) + "\n") + PY + - name: Publish complete prediction seal before target scoring + if: env.EVALUATION_MODE == 'evaluate' + uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.0 + with: + name: tracking-cloth-prediction-seal-${{ github.run_id }}-${{ github.run_attempt }} + path: | + ${{ env.CLOTH_OUT }}/protocol.json + ${{ env.CLOTH_OUT }}/source_fit.json + ${{ env.CLOTH_OUT }}/dataset_manifest.json + ${{ env.CLOTH_OUT }}/prediction_seal.json + ${{ env.CLOTH_OUT }}/run_manifest.json + ${{ env.CLOTH_OUT }}/evaluation_request.json + ${{ env.CLOTH_OUT }}/DATA_LICENSE.txt + retention-days: 90 + if-no-files-found: error + - name: Score only sealed twisting forecasts + if: env.EVALUATION_MODE == 'evaluate' + shell: bash + run: | + set -euo pipefail + "$CLOTH_PY" -m experiments.tracking_cloth_deformation_v1.run \ + --dataset-root "$DATASET_ROOT" --output "$CLOTH_OUT" --stage score + - name: Publish operator summary + if: always() && env.CLOTH_OUT != '' + shell: bash + run: | + if test -f "$CLOTH_OUT/report.md"; then cat "$CLOTH_OUT/report.md" >> "$GITHUB_STEP_SUMMARY"; fi + if test -f "$CLOTH_OUT/failure.json"; then + echo '## Incomplete run: no scientific conclusion' >> "$GITHUB_STEP_SUMMARY" + cat "$CLOTH_OUT/failure.json" >> "$GITHUB_STEP_SUMMARY" + fi + - name: Upload aggregate evidence only, never raw recordings or trajectory arrays + if: always() && env.CLOTH_OUT != '' + uses: actions/upload-artifact@043fb46d1a93c77aae656e7c1c64a875d1fc6a0a # v7.0.0 + with: + name: tracking-cloth-${{ env.EVALUATION_MODE }}-${{ github.run_id }}-${{ github.run_attempt }} + path: | + ${{ env.CLOTH_OUT }}/*.json + ${{ env.CLOTH_OUT }}/*.csv + ${{ env.CLOTH_OUT }}/report.md + ${{ env.CLOTH_OUT }}/DATA_LICENSE.txt + retention-days: 90 + if-no-files-found: warn + - name: Remove private predictions and isolated environment + if: always() + shell: bash + run: | + if test -n "${CLOTH_OUT:-}"; then rm -rf "$CLOTH_OUT/private_predictions"; fi + if test -n "${CLOTH_VENV:-}"; then rm -rf "$CLOTH_VENV"; fi diff --git a/docs/workflow_inventory_budget.md b/docs/workflow_inventory_budget.md index 2e092f9ed..ba119f4a9 100644 --- a/docs/workflow_inventory_budget.md +++ b/docs/workflow_inventory_budget.md @@ -26,7 +26,7 @@ At source revision The machine-readable contract is `.github/quality/workflow-inventory-budget-v1.json`. -## Completed one-shot retirement and current reviewed increase +## Completed one-shot retirement and reviewed permanent increases The twelve historical one-shot files are absent from `.github/workflows`. Their exact Git blobs are retained below @@ -46,18 +46,27 @@ ordinary workflows with zero temporary-looking files. The permanent addition: it runs the target-free controlled falsification study, requires primary/replay byte identity, binds regenerated output to the retained result, and records the exact reviewed head and canonical Python/NumPy runtime. Its -addition raises the checked-in ceiling by one, from 81 to 82, without changing +addition raised the checked-in ceiling by one, from 81 to 82, without changing the retirement target. +The permanent `tracking-cloth-evaluation.yml` workflow is a second deliberately +reviewed addition. It provides synthetic contract validation plus a read-only, +main-branch-only public real-data pilot on the protected `gpuserver4090` runner. +The workflow freezes source-only shake fitting before twisting-target scoring, +publishes the complete prediction seal before scoring, retains aggregate +evidence only, and requires a canonical single-file authorization commit for +the first target-scored run. Its addition raises the checked-in ceiling by one, +from 82 to 83, without changing the retirement target. + The exact active inventory and targets are therefore: -- 82 checked-in workflows; +- 83 checked-in workflows; - zero temporary-looking workflow files; - a retirement target of at most 81 checked-in workflows; and - a retirement target of zero temporary-looking workflows. -The one-workflow retirement gap is intentional and visible. A future -consolidation or retirement should lower the ceiling back to 81 in the same +The two-workflow retirement gap is intentional and visible. Future +consolidations or retirements should lower the ceiling toward 81 in the same change rather than silently reusing that capacity. Validate the active inventory with: diff --git a/experiments/tracking_cloth_deformation_v1/README.md b/experiments/tracking_cloth_deformation_v1/README.md new file mode 100644 index 000000000..479e138dd --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/README.md @@ -0,0 +1,182 @@ +# Tracking Cloth Deformation: shake-to-twist evaluation pilot + +A maintained workflow for the user-installed **real public motion-capture data** +from Coltraro, Borras, Alberich-Carraminana and Torras (2025), +*Tracking cloth deformation: A novel dataset for closing the sim-to-real gap for +robotic cloth manipulation learning*, DOI `10.1177/02783649251317617`. +Dataset: . + +This is an executable, explicitly limited **physical-baseline qualification and +public-data pilot**, not a reproduction of PhysTwin or clothilde-sim, not a +paper-ready validation of the complete BayesianPhysTwin API, and not another +request for new physical recordings. It changes no historical result or claim. + +## Run on GitHub + +Open **Actions -> Tracking Cloth Deformation evaluation -> Run workflow**. +The evaluation job requires labels `[self-hosted, Linux, X64, gpuserver6000]`. +The intended server is workstation2; routing follows the user-specified label, +not the earlier gpuserver4090 installation note. The default read-only cache is: + +```text +/home/github-runner/.cache/datasets/tracking-cloth-deformation-v1-zenodo-14644526 +``` + +The modes are: + +| Mode | Data access | Outputs | +| --- | --- | --- | +| `inventory` | Archive integrity/hashes and file names; no numeric trajectories | Verified roster, licensing record, provenance | +| `source_only` (default) | 32 shaking recordings only | Four-fold source scores, parameter weights, empirical guards, variance calibration | +| `evaluate` | Source fitting, then target initialization/corners; all 32 predictions sealed before target free-marker scoring | Per-record/per-specimen scores, paired contrasts, coverage/width, fallback audit | + +Choose `source_only` first. Review source model competence and the initialization +geometry/units assumptions before explicitly choosing `evaluate`. The latter +repeats the same source fit from the same versioned protocol, publishes its +prediction-seal artifact, and only then runs a separate target-scoring process. +A failed preparation or upload prevents scoring. A changed source model, +protocol, implementation, dataset or prediction artifact invalidates the seal. +All runs remain pilots: prior public outcome exposure is unknown, and rerunning +already-scored targets never creates fresh confirmation. + +The cache is never downloaded, extracted, modified, or chmod'ed. A private +scratch directory and isolated Python environment are created in `RUNNER_TEMP`. +Only NumPy and CPU are required. Four CPU workers are the default. Pull requests +run **synthetic tests on GitHub-hosted runners only**; they never access the +self-hosted dataset. No runner/organization secret or repository-write permission +is needed. Concurrency does not cancel an active evidence run. + +## Frozen v1 design + +`protocol.json` fixes the following before target numerical use: + +* Exact Zenodo record and published archive MD5 + `b4868b702f8a42b2ea1069d0f1a3b8f6`; 120 extracted CSVs must byte-match the ZIP. +* The complete free-hanging factorial: 4 materials x 2 sizes x 2 speeds x + 2 grasp modes x 2 motions. All 32 shaking recordings are sources and all 32 + twisting recordings are targets. The other 56 collision recordings remain + unused, not opportunistic replacements. +* A 1-second all-marker initialization prefix and a 5-second scored forecast. + Data are regularly subsampled from 120 Hz to 30 Hz; dynamics use eight + integration substeps. No time alignment is fitted to target outcomes. +* The future measured trajectories of the **two driven corners are supplied as + prescribed boundary conditions**, as in the dataset's simulator use case. + They are not logged robot commands and are not scored. No future free-marker + coordinate enters fitting, state initialization, model selection or prediction. +* Four-fold leave-one-speed/grasp-recording-out source fitting within each + material-size specimen. A candidate guard is fixed from these source folds, + not selected from twisting scores. + +CSV parsing follows the paper's numeric `Frame, Time, X1,Y1,Z1,...` layout and +preserves missing-marker masks. Only source initialization geometry is used to +resolve metre/centimetre/millimetre scale; all source units must agree, and that +scale is frozen for targets. The pilot assumes an initially near-vertical regular +marker mesh: A2 has five rows of four markers, A3 four rows of three. Geometry +is ordered from the initial frame, with the upper-row endpoints as driven +corners. This is **not asserted to reproduce the supplied MATLAB ordering**. +Unsupported or ambiguous initialization fails the run, with no dropped cases or +outcome-dependent remapping. The four supplied MATLAB readers are hashed but not +executed. Prefix and corner gaps use causal carry-forward, never interpolation +from future free-marker observations. Missing free-marker ground truth is omitted +with the same mask for every arm and its count is reported. + +## What the pilot actually computes + +The physical backend is a small equal-marker-mass 3-D spring mesh with structural, +shear and two-hop bending springs, gravity, viscous damping, symplectic +integration and prescribed corner positions. Its rest lengths come from the +initial frame. The fixed bank contains stiffness-per-mass values +`[100, 400, 1600]` and damping-per-mass values `[0.5, 2, 8]`. This is a transparent +qualification baseline, not a high-fidelity material/contact model. + +The seven arms are: + +| Arm | Definition | +| --- | --- | +| `persistence` | Hold the last permitted prefix free-marker position | +| `nominal_physics` | Nominal `(400,2)` spring model from initial state, no prefix-end reset | +| `last_residual` | Same nominal forecast plus its last prefix residual, outside simulator state | +| `nominal_state_injection` | Nominal model reset to prefix-end positions and backward-only velocity | +| `map_physics` | Source-selected maximum-weight parameter member with the same reset | +| `bayesian_physics` | Source-weighted parameter-bank average with the same reset | +| `guarded_bayesian_physics` | Complete candidate mean/variance or exact nominal mean/variance | + +The source model weights are an explicit **generalized Bayesian/Gibbs** update: +`w(k) proportional to exp(-sum_record MSE(record,k)/(2*temperature))`. +Each source recording contributes one normalized loss, not thousands of +independent marker likelihoods. Temperature is the source-only median best-member +MSE, floored at `1 mm` squared. These are not claimed to be calibrated physical +parameter posterior probabilities. Parameter uncertainty contributes the +between-member variance; three source-only out-of-fold horizon bins supply a +residual variance floor. Scoring uses a diagonal moment-matched Gaussian, not a +claim of calibrated joint trajectory covariance. + +The empirical source guard accepts a specimen only when its out-of-fold +candidate mean RMSE beats **both** nominal physics and last-residual by at least +1%, and none of its four source folds regresses versus nominal physics. This is +a conservative source selection rule, **not a finite-sample deployment safety +certificate**. The unguarded controls are always reported. A rejected candidate +reuses the exact nominal mean and variance, and scoring verifies equality. + +## Results and interpretation + +Primary endpoint: free-marker Euclidean trajectory RMSE in millimetres, averaged +within each recording, then equally over its four speed/grasp conditions and +then over eight material-size specimens. This is **not** the dataset paper's +mass-matrix norm; no unverified marker mass quadrature is invented. + +Supporting endpoints are mean marker Euclidean error, coordinate Gaussian NLL, +90% marginal coverage **with interval width**, missing scored-marker counts, +accepted/rejected recordings, harmful accepted records, worst-specimen regret, +and exact fallback violations. Comparisons include the strong last-residual and +MAP controls, not just the nominal simulator. + +Paired bootstrap intervals resample the eight specimens. A four-material +cluster sensitivity analysis is also reported because sizes may not constitute +independent draws of material properties. These are small-sample, exploratory, +non-simultaneous intervals, not a multi-comparison confirmation certificate. +Frames, marker coordinates and the 32 recordings are not treated as independent +physical specimens. No case is silently removed from the registered factorial. + +Artifacts contain `protocol.json`, dataset/code hashes, `source_fit.json`, +`source_scores.csv`, `prediction_seal.json`, `target_scores.csv`, +`specimen_scores.csv`, `metrics.json`, `run_manifest.json`, `report.md`, and the +included license. Prediction trajectory arrays stay in private runner scratch +and are **not uploaded**, as are all raw recordings and the ZIP. Reports are +retained for 90 days, not forever; a paper claim needs a separate durable evidence +intake after scientific review. A technical failure emits `failure.json` and no +complete scientific decision. Green unit tests are synthetic software evidence. + +## License conflict + +The user-retained Zenodo metadata says CC BY-SA 4.0, whereas the included +`License.txt` says **CC BY-NC-SA 4.0**. The stricter included noncommercial policy +is retained pending author clarification. The archive, metadata and originals +remain untouched in the cache; the exact license and hashes accompany output. +This workflow does not resolve the conflict or grant commercial/redistribution +rights. Do not relicense the dataset as the repository's code license. + +## Local equivalent + +From the repository root, in an environment with the pinned requirements: + +```bash +python -m experiments.tracking_cloth_deformation_v1.run \ + --dataset-root /home/github-runner/.cache/datasets/tracking-cloth-deformation-v1-zenodo-14644526 \ + --output /tmp/tracking-cloth-source-v1 --stage source --workers 4 +``` + +For an explicitly authorized pilot, use `--stage predict` in a **fresh output +directory**, preserve the complete prediction seal, and then run `--stage score` +against that same directory. The workflow additionally uploads the seal before +scoring. Commands never revise `claims.json` or manuscript prose. + +## Workflow lifecycle and evidence classification + +Classification: **operational prerequisite and scored diagnostic** for the +untested external real-cloth shake-to-twist question, specifically authorized by +the user after installing this dataset. The existing Cloth Sim2Real workflow is +bound to a different Zenodo release, simulator, cohort and historical evidence; +repurposing it would silently change its scientific contract. This one maintained +entry point provides inventory, source qualification and sealed scoring, with +experiment logic in versioned Python. No one-shot bootstrap workflow is added. diff --git a/experiments/tracking_cloth_deformation_v1/__init__.py b/experiments/tracking_cloth_deformation_v1/__init__.py new file mode 100644 index 000000000..209c3cfcc --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/__init__.py @@ -0,0 +1 @@ +"""Private, public-data cloth forecasting pilot (not a supported package API).""" diff --git a/experiments/tracking_cloth_deformation_v1/active_probe.py b/experiments/tracking_cloth_deformation_v1/active_probe.py new file mode 100644 index 000000000..bb641789f --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/active_probe.py @@ -0,0 +1,351 @@ +"""Outcome-blind finite-model probe selection for the tracking-cloth pilot. + +The selection rules consume only source-frozen model-disagreement matrices and +current discrete-model weights. Outcomes are requested from the supplied +mapping only after an action has been selected. This is a retrospective replay +of already recorded actions, not an online robot controller or a safety policy. +""" + +from __future__ import annotations + +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from typing import Any + +import numpy as np +from numpy.typing import NDArray + +FloatArray = NDArray[np.float64] + +POLICIES = ("fixed_order", "parameter_information", "task_directed") + + +def _finite_array(value: object, *, name: str, ndim: int | None = None) -> FloatArray: + array = np.asarray(value, dtype=np.float64).copy() + if ndim is not None and array.ndim != ndim: + raise ValueError(f"{name} must have {ndim} dimensions") + if not np.all(np.isfinite(array)): + raise ValueError(f"{name} must contain only finite values") + return array + + +def normalize_weights(value: object) -> FloatArray: + """Return a finite, strictly positive normalized model-weight vector.""" + weights = _finite_array(value, name="weights", ndim=1) + if weights.size < 2 or np.any(weights < 0.0) or not np.any(weights > 0.0): + raise ValueError("weights must be a nonnegative vector with positive mass") + total = float(weights.sum()) + if not np.isfinite(total) or total <= 0.0: + raise ValueError("weights must have finite positive mass") + weights /= total + tiny = np.finfo(np.float64).tiny + weights = np.maximum(weights, tiny) + weights /= weights.sum() + weights.setflags(write=False) + return weights + + +def _temperature(value: float) -> float: + result = float(value) + if not np.isfinite(result) or result <= 0.0: + raise ValueError("temperature must be finite and positive") + return result + + +def update_weights(weights: object, loss: object, temperature: float) -> FloatArray: + """Apply one Gibbs update from one normalized trajectory loss vector.""" + prior = normalize_weights(weights) + losses = _finite_array(loss, name="loss", ndim=1) + if losses.shape != prior.shape or np.any(losses < 0.0): + raise ValueError("loss must be nonnegative and match weights") + temperature = _temperature(temperature) + logits = np.log(prior) - losses / (2.0 * temperature) + logits -= float(np.max(logits)) + return normalize_weights(np.exp(logits)) + + +def weights_from_records(losses: object, temperature: float) -> FloatArray: + """Fit one equal-record generalized-Bayes prior over model members.""" + matrix = _finite_array(losses, name="losses", ndim=2) + if matrix.shape[0] < 1 or matrix.shape[1] < 2 or np.any(matrix < 0.0): + raise ValueError( + "losses must have shape (records, models), with nonnegative values" + ) + temperature = _temperature(temperature) + return update_weights(np.ones(matrix.shape[1]), matrix.sum(axis=0), temperature) + + +def pairwise_trajectory_mse(bank: object, valid: object) -> FloatArray: + """Pairwise mean squared Euclidean trajectory disagreement. + + ``bank`` has shape ``(models, time, points, 3)`` and ``valid`` has shape + ``(time, points)``. The metric matches the point-wise squared-error scale + used by the spring-pilot likelihood. + """ + values = _finite_array(bank, name="bank", ndim=4) + if values.shape[0] < 2 or values.shape[-1] != 3: + raise ValueError("bank must have shape (models>=2, time, points, 3)") + mask = np.asarray(valid, dtype=bool) + if mask.shape != values.shape[1:3] or not np.any(mask): + raise ValueError("valid must select at least one time-point entry") + selected = values[:, mask, :] + delta = selected[:, None, :, :] - selected[None, :, :, :] + distance = np.mean(np.sum(delta * delta, axis=-1), axis=-1) + distance = 0.5 * (distance + distance.T) + np.fill_diagonal(distance, 0.0) + if np.any(distance < -1e-15): + raise ValueError("pairwise disagreement must be nonnegative") + distance = np.maximum(distance, 0.0) + distance.setflags(write=False) + return distance + + +def validate_distance_matrix(value: object, *, models: int, name: str) -> FloatArray: + matrix = _finite_array(value, name=name, ndim=2) + if matrix.shape != (models, models): + raise ValueError(f"{name} must have shape ({models}, {models})") + scale = max(float(np.max(np.abs(matrix))), 1.0) + if not np.allclose(matrix, matrix.T, atol=1e-12 * scale, rtol=1e-10): + raise ValueError(f"{name} must be symmetric") + if not np.allclose(np.diag(matrix), 0.0, atol=1e-12 * scale, rtol=0.0): + raise ValueError(f"{name} must have a zero diagonal") + if np.any(matrix < -1e-12 * scale): + raise ValueError(f"{name} must be nonnegative") + matrix = np.maximum(0.5 * (matrix + matrix.T), 0.0) + matrix.setflags(write=False) + return matrix + + +def entropy(weights: object) -> float: + probabilities = normalize_weights(weights) + return float(-np.sum(probabilities * np.log(probabilities))) + + +def model_spread(weights: object, distance: object) -> float: + """Expected squared pairwise disagreement divided by two.""" + probabilities = normalize_weights(weights) + matrix = validate_distance_matrix( + distance, models=probabilities.size, name="distance" + ) + return float(0.5 * probabilities @ matrix @ probabilities) + + +def pseudo_posteriors( + weights: object, probe_distance: object, temperature: float +) -> FloatArray: + """Posterior for each model-member pseudo-outcome under source noise. + + Row ``j`` is the posterior obtained if member ``j`` generated the noiseless + probe mean. The source-frozen temperature softens pairwise separation. + """ + probabilities = normalize_weights(weights) + distance = validate_distance_matrix( + probe_distance, models=probabilities.size, name="probe_distance" + ) + temperature = _temperature(temperature) + rows = np.stack( + [ + update_weights(probabilities, distance[index], temperature) + for index in range(probabilities.size) + ] + ) + rows.setflags(write=False) + return rows + + +def parameter_information_utility( + weights: object, probe_distance: object, temperature: float +) -> float: + """Expected entropy reduction over discrete physical-model members.""" + probabilities = normalize_weights(weights) + posteriors = pseudo_posteriors(probabilities, probe_distance, temperature) + expected = sum( + float(probabilities[index]) * entropy(posteriors[index]) + for index in range(probabilities.size) + ) + return float(max(entropy(probabilities) - expected, 0.0)) + + +def task_variance_reduction_utility( + weights: object, + probe_distance: object, + target_distance: object, + temperature: float, +) -> float: + """Expected fractional contraction of held-out-task model spread.""" + probabilities = normalize_weights(weights) + target = validate_distance_matrix( + target_distance, models=probabilities.size, name="target_distance" + ) + current = model_spread(probabilities, target) + if current <= np.finfo(np.float64).eps: + return 0.0 + posteriors = pseudo_posteriors(probabilities, probe_distance, temperature) + expected = sum( + float(probabilities[index]) * model_spread(posteriors[index], target) + for index in range(probabilities.size) + ) + reduction = (current - expected) / current + return float(np.clip(reduction, 0.0, 1.0)) + + +def select_action( + *, + policy: str, + weights: object, + remaining_actions: Sequence[str], + probe_distances: Mapping[str, object], + target_distance: object, + temperature: float, + fixed_order: Sequence[str], +) -> tuple[str, dict[str, float | None]]: + """Select one action without consuming an action outcome.""" + if policy not in POLICIES: + raise ValueError(f"unknown policy: {policy}") + remaining = tuple(str(action) for action in remaining_actions) + if not remaining or len(remaining) != len(set(remaining)): + raise ValueError("remaining_actions must be nonempty and unique") + if set(probe_distances) != set(fixed_order): + raise ValueError( + "probe_distances and fixed_order must define the same action roster" + ) + if not set(remaining).issubset(probe_distances): + raise ValueError("remaining action is absent from probe_distances") + probabilities = normalize_weights(weights) + target = validate_distance_matrix( + target_distance, models=probabilities.size, name="target_distance" + ) + utilities: dict[str, float | None] = {} + if policy == "fixed_order": + order = tuple(str(action) for action in fixed_order) + if len(order) != len(set(order)) or set(order) != set(probe_distances): + raise ValueError("fixed_order must contain every action exactly once") + selected = next(action for action in order if action in remaining) + utilities.update({action: None for action in remaining}) + return selected, utilities + for action in sorted(remaining): + probe = validate_distance_matrix( + probe_distances[action], + models=probabilities.size, + name=f"probe_distances[{action}]", + ) + if policy == "parameter_information": + utility = parameter_information_utility(probabilities, probe, temperature) + else: + utility = task_variance_reduction_utility( + probabilities, probe, target, temperature + ) + utilities[action] = utility + maximum = max(float(value) for value in utilities.values() if value is not None) + tolerance = 1e-14 * max(1.0, abs(maximum)) + selected = min( + action + for action, value in utilities.items() + if value is not None and maximum - float(value) <= tolerance + ) + return selected, utilities + + +@dataclass(frozen=True) +class PolicyState: + budget: int + selected_actions: tuple[str, ...] + weights: FloatArray + steps: tuple[dict[str, Any], ...] + + def __post_init__(self) -> None: + if ( + isinstance(self.budget, bool) + or int(self.budget) != self.budget + or self.budget < 0 + ): + raise ValueError("budget must be a nonnegative integer") + weights = normalize_weights(self.weights) + actions = tuple(str(action) for action in self.selected_actions) + if len(actions) != len(set(actions)) or len(actions) != self.budget: + raise ValueError("selected_actions must be unique and match budget") + object.__setattr__(self, "budget", int(self.budget)) + object.__setattr__(self, "selected_actions", actions) + object.__setattr__(self, "weights", weights) + object.__setattr__(self, "steps", tuple(dict(step) for step in self.steps)) + + +def simulate_policy( + *, + policy: str, + initial_weights: object, + probe_distances: Mapping[str, object], + target_distance: object, + observed_losses: Mapping[str, object], + temperature: float, + fixed_order: Sequence[str], + budgets: Sequence[int], +) -> dict[int, PolicyState]: + """Replay a sequential policy and consume only selected probe outcomes.""" + requested = tuple(int(value) for value in budgets) + action_count = len(probe_distances) + if ( + not requested + or tuple(sorted(set(requested))) != requested + or requested[0] != 0 + or requested[-1] > action_count + ): + raise ValueError("budgets must be sorted unique values starting at zero") + if set(observed_losses) != set(probe_distances): + raise ValueError("observed_losses must expose the complete registered roster") + weights = normalize_weights(initial_weights) + states: dict[int, PolicyState] = {0: PolicyState(0, (), weights, ())} + selected: list[str] = [] + steps: list[dict[str, Any]] = [] + remaining = set(probe_distances) + for step_index in range(1, requested[-1] + 1): + action, utilities = select_action( + policy=policy, + weights=weights, + remaining_actions=tuple(sorted(remaining)), + probe_distances=probe_distances, + target_distance=target_distance, + temperature=temperature, + fixed_order=fixed_order, + ) + entropy_before = entropy(weights) + target_before = model_spread(weights, target_distance) + # This is the only point at which the selected action's outcome is read. + selected_loss = observed_losses[action] + weights = update_weights(weights, selected_loss, temperature) + selected.append(action) + remaining.remove(action) + steps.append( + { + "step": step_index, + "selected_action": action, + "utilities": utilities, + "entropy_before": entropy_before, + "entropy_after": entropy(weights), + "target_model_spread_before": target_before, + "target_model_spread_after": model_spread(weights, target_distance), + } + ) + if step_index in requested: + states[step_index] = PolicyState( + step_index, tuple(selected), weights, tuple(steps) + ) + return states + + +__all__ = [ + "POLICIES", + "PolicyState", + "entropy", + "model_spread", + "normalize_weights", + "pairwise_trajectory_mse", + "parameter_information_utility", + "pseudo_posteriors", + "select_action", + "simulate_policy", + "task_variance_reduction_utility", + "update_weights", + "validate_distance_matrix", + "weights_from_records", +] diff --git a/experiments/tracking_cloth_deformation_v1/active_probe_protocol.json b/experiments/tracking_cloth_deformation_v1/active_probe_protocol.json new file mode 100644 index 000000000..8ee3f9b47 --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/active_probe_protocol.json @@ -0,0 +1,53 @@ +{ + "study_id": "tracking-cloth-task-directed-active-probe-v1", + "dataset_record": "14644526", + "archive_md5": "b4868b702f8a42b2ea1069d0f1a3b8f6", + "csv_count": 120, + "source_motion": "shake", + "target_motion": "twist", + "materials": ["cotton", "denim", "polyester", "wool"], + "sizes": ["A2", "A3"], + "speeds": ["fast", "slow"], + "grasps": ["hands", "hanger"], + "sample_stride": 4, + "prefix_seconds": 1.0, + "forecast_seconds": 5.0, + "initial_complete_frame_deadline_seconds": 0.25, + "stiffness_per_mass": [100.0, 400.0, 1600.0], + "damping_per_mass": [0.5, 2.0, 8.0], + "nominal_parameters": [400.0, 2.0], + "integration_substeps": 8, + "gravity_m_s2": 9.80665, + "measurement_floor_m": 0.001, + "fold_rule": "leave-one-material-out", + "probe_conditions": [ + "fast_hands", + "fast_hanger", + "slow_hands", + "slow_hanger" + ], + "fixed_probe_order": [ + "fast_hands", + "fast_hanger", + "slow_hands", + "slow_hanger" + ], + "probe_policies": [ + "fixed_order", + "parameter_information", + "task_directed" + ], + "probe_budgets": [0, 1, 2, 4], + "primary_budget": 1, + "selection_model": "finite-model pseudo-outcome expected contraction", + "selection_templates": "other-material shake and twist-input model disagreement only", + "held_material_candidate_inputs_used_for_selection": false, + "held_material_twist_inputs_used_for_selection": false, + "bootstrap_repetitions": 10000, + "bootstrap_seed": 14644526, + "information_regime": "recorded-corner-boundary-conditioned target forecast; source-frozen probe selection", + "evidence_class": "cross-validated public-real-data active-probe pilot; prior exposure not assumed absent", + "license_policy": "CC-BY-NC-SA-4.0 pending author clarification of conflicting CC-BY-SA-4.0 metadata", + "raw_data_upload": false, + "paper_claim_authorized": false +} diff --git a/experiments/tracking_cloth_deformation_v1/active_probe_run.py b/experiments/tracking_cloth_deformation_v1/active_probe_run.py new file mode 100644 index 000000000..ab0bc3c23 --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/active_probe_run.py @@ -0,0 +1,344 @@ +"""Belief construction and sealing helpers for the active-probe cloth pilot. + +The helpers in this module are deliberately outcome-blind. They turn frozen +finite-model weights into complete mean/variance trajectories and expose the +small set of deterministic operations needed by the registered retrospective +probe study. They do not select from or inspect held-out twisting outcomes. +""" + +from __future__ import annotations + +import hashlib +import json +from collections.abc import Mapping, Sequence +from typing import Any + +import numpy as np + +from .active_probe import POLICIES, normalize_weights +from .data import Inputs +from .model import Predictions, horizon_bins + + +def _finite_array(value: object, *, name: str, ndim: int | None = None) -> np.ndarray: + array = np.asarray(value, dtype=np.float64).copy() + if ndim is not None and array.ndim != ndim: + raise ValueError(f"{name} must have {ndim} dimensions") + if not np.all(np.isfinite(array)): + raise ValueError(f"{name} must contain only finite values") + return array + + +def _residual_variance(value: object, *, name: str) -> np.ndarray: + variance = _finite_array(value, name=name, ndim=1) + if variance.shape != (3,) or np.any(variance <= 0.0): + raise ValueError(f"{name} must contain three positive horizon variances") + return variance + + +def _variance_field( + inputs: Inputs, shape: tuple[int, ...], residual_variance: object +) -> np.ndarray: + residual = _residual_variance(residual_variance, name="residual_variance") + bins = horizon_bins(inputs) + variance = np.broadcast_to(residual[bins, None, None], shape).copy() + if np.any(variance <= 0.0) or not np.all(np.isfinite(variance)): + raise ValueError("predictive variance must be finite and positive") + return variance + + +def active_mask(inputs: Inputs) -> np.ndarray: + """Return the registered post-prefix, non-driven-marker mask.""" + times = np.asarray(inputs.times) + order = np.asarray(inputs.order) + corners = np.asarray(inputs.corners, dtype=int) + if times.ndim != 1 or order.ndim != 1: + raise ValueError("times and marker order must be one-dimensional") + if not 0 <= int(inputs.cutoff) < len(times): + raise ValueError("cutoff is outside the trajectory") + if corners.shape != (2,) or np.any(corners < 0) or np.any(corners >= len(order)): + raise ValueError("exactly two valid driven-corner indices are required") + mask = np.ones((len(times), len(order)), dtype=bool) + mask[: int(inputs.cutoff) + 1] = False + mask[:, corners] = False + if not np.any(mask): + raise ValueError("active mask contains no scored samples") + return mask + + +def array_digest(value: object) -> str: + """Hash array values together with their canonical shape and dtype.""" + array = np.ascontiguousarray(np.asarray(value)) + metadata = json.dumps( + {"dtype": array.dtype.str, "shape": list(array.shape)}, + sort_keys=True, + separators=(",", ":"), + ).encode("utf-8") + digest = hashlib.sha256() + digest.update(metadata) + digest.update(b"\0") + digest.update(array.tobytes(order="C")) + return digest.hexdigest() + + +def belief_digest(mean: object, variance: object) -> str: + """Bind a complete mean and variance trajectory into one content ID.""" + payload = json.dumps( + {"mean": array_digest(mean), "variance": array_digest(variance)}, + sort_keys=True, + separators=(",", ":"), + ).encode("utf-8") + return hashlib.sha256(payload).hexdigest() + + +def weighted_belief( + prediction: Predictions, weights: object, residual_variance: object +) -> tuple[np.ndarray, np.ndarray]: + """Moment-match one finite-model trajectory mixture.""" + bank = _finite_array(prediction.bank, name="prediction.bank", ndim=4) + if bank.shape[0] < 2 or bank.shape[-1] != 3: + raise ValueError("prediction bank must have shape (models, time, markers, 3)") + probabilities = normalize_weights(weights) + if probabilities.shape != (bank.shape[0],): + raise ValueError("weights must match the prediction-bank model count") + mean = np.einsum("k,ktnd->tnd", probabilities, bank) + variance = np.einsum("k,ktnd->tnd", probabilities, (bank - mean[None, ...]) ** 2) + variance += _variance_field(prediction.inputs, mean.shape, residual_variance) + if not np.all(np.isfinite(mean)) or np.any(variance <= 0.0): + raise ValueError("invalid moment-matched belief") + return mean, variance + + +def loss_vector(prediction: Predictions, truth: object) -> np.ndarray: + """Return one normalized trajectory loss for every model member.""" + bank = _finite_array(prediction.bank, name="prediction.bank", ndim=4) + observed = np.asarray(truth, dtype=np.float64) + if observed.shape != bank.shape[1:]: + raise ValueError("truth must match one prediction-bank trajectory") + valid = active_mask(prediction.inputs) & np.isfinite(observed).all(axis=2) + if not np.any(valid): + raise ValueError("truth contains no evaluable active samples") + if not np.all(np.isfinite(bank[:, valid, :])): + raise ValueError("prediction bank is nonfinite on evaluable samples") + error = bank[:, valid, :] - observed[valid][None, ...] + losses = np.mean(np.sum(error * error, axis=-1), axis=1) + losses = np.maximum(losses, 0.0) + losses.setflags(write=False) + return losses + + +def posterior_temperature(losses: object, measurement_floor_m: float) -> float: + """Freeze the generalized-Bayes temperature from source records only.""" + matrix = _finite_array(losses, name="losses", ndim=2) + if matrix.shape[0] < 1 or matrix.shape[1] < 2 or np.any(matrix < 0.0): + raise ValueError("losses must be nonnegative with shape (records, models>=2)") + floor = float(measurement_floor_m) + if not np.isfinite(floor) or floor <= 0.0: + raise ValueError("measurement_floor_m must be finite and positive") + return max(float(np.median(np.min(matrix, axis=1))), floor * floor) + + +def _last_residual_mean( + prediction: Predictions, prefix_last: object, boundary: object +) -> np.ndarray: + prefix = _finite_array(prefix_last, name="prefix_last", ndim=2) + if prefix.shape != prediction.nominal.shape[1:]: + raise ValueError("prefix_last must match the marker-state shape") + driven = _finite_array(boundary, name="boundary", ndim=3) + if driven.shape != (prediction.nominal.shape[0], 2, 3): + raise ValueError("boundary must have shape (time, 2, 3)") + mean = _finite_array(prediction.nominal, name="prediction.nominal", ndim=3) + residual = prefix - mean[int(prediction.inputs.cutoff)] + mean += residual[None, ...] + mean[:, np.asarray(prediction.inputs.corners, dtype=int)] = driven + return mean + + +def calibrated_residuals( + records: Sequence[tuple[Predictions, object]], + weights: object, + protocol: Mapping[str, Any], +) -> dict[str, list[float]]: + """Estimate three equal-record residual-variance bins from source data.""" + if not records: + raise ValueError("at least one source record is required") + floor = float(protocol["measurement_floor_m"]) + if not np.isfinite(floor) or floor <= 0.0: + raise ValueError("measurement_floor_m must be finite and positive") + floor2 = floor * floor + pooled: dict[str, list[list[float]]] = { + "bayesian": [[], [], []], + "nominal_physics": [[], [], []], + "last_residual": [[], [], []], + } + for prediction, raw_truth in records: + truth = np.asarray(raw_truth, dtype=np.float64) + if truth.shape != prediction.nominal.shape: + raise ValueError("source truth must match the prediction trajectory") + valid = active_mask(prediction.inputs) & np.isfinite(truth).all(axis=2) + if not np.any(valid): + raise ValueError("source record contains no evaluable samples") + probabilities = normalize_weights(weights) + if probabilities.shape != (prediction.bank.shape[0],): + raise ValueError("weights must match every source prediction bank") + bayesian_mean = np.einsum("k,ktnd->tnd", probabilities, prediction.bank) + ensemble_variance = np.einsum( + "k,ktnd->tnd", + probabilities, + (prediction.bank - bayesian_mean[None, ...]) ** 2, + ) + last_residual = _last_residual_mean( + prediction, + prediction.inputs.prefix[-1], + prediction.inputs.boundary, + ) + means = { + "bayesian": bayesian_mean, + "nominal_physics": prediction.nominal, + "last_residual": last_residual, + } + bins = horizon_bins(prediction.inputs) + for name, mean in means.items(): + error2 = (mean - truth) ** 2 + if name == "bayesian": + error2 = error2 - ensemble_variance + for bin_index in range(3): + selected = valid & (bins[:, None] == bin_index) + if not np.any(selected): + raise ValueError("empty source calibration horizon bin") + pooled[name][bin_index].append( + max(float(np.mean(error2[selected])), floor2) + ) + return { + name: [max(float(np.mean(values)), floor2) for values in horizon_values] + for name, horizon_values in pooled.items() + } + + +def validate_protocol(protocol: Mapping[str, Any]) -> None: + """Validate the active-probe roster and its information boundary.""" + conditions = tuple(str(value) for value in protocol["probe_conditions"]) + fixed_order = tuple(str(value) for value in protocol["fixed_probe_order"]) + policies = tuple(str(value) for value in protocol["probe_policies"]) + budgets = tuple(protocol["probe_budgets"]) + if len(conditions) != 4 or len(set(conditions)) != len(conditions): + raise ValueError("probe_conditions must contain four unique actions") + if fixed_order != conditions: + raise ValueError("fixed_probe_order must equal the registered action roster") + if policies != POLICIES: + raise ValueError("probe_policies must match the implemented policy roster") + if budgets != (0, 1, 2, 4): + raise ValueError("probe_budgets must be exactly [0, 1, 2, 4]") + if protocol["primary_budget"] not in budgets: + raise ValueError("primary_budget must be registered") + if protocol["held_material_candidate_inputs_used_for_selection"] is not False: + raise ValueError("held-material candidate inputs may not drive selection") + if protocol["held_material_twist_inputs_used_for_selection"] is not False: + raise ValueError("held-material twist inputs may not drive selection") + if protocol["paper_claim_authorized"] is not False: + raise ValueError("the public-data pilot may not self-authorize a paper claim") + + +def arm_specs(protocol: Mapping[str, Any]) -> tuple[str, ...]: + """Return the complete 18-arm registered comparison roster.""" + validate_protocol(protocol) + controls = ("nominal_physics", "last_residual") + single_probes = tuple( + f"single_probe_{condition}" for condition in protocol["probe_conditions"] + ) + policies = tuple( + f"{policy}_k{budget}" + for policy in protocol["probe_policies"] + for budget in protocol["probe_budgets"] + ) + result = controls + single_probes + policies + if len(result) != 18 or len(set(result)) != len(result): + raise ValueError("active-probe arm roster is not complete and unique") + return result + + +def _mapping(value: object, *, name: str) -> Mapping[str, Any]: + if not isinstance(value, Mapping): + raise ValueError(f"{name} must be a mapping") + return value + + +def build_belief_arms( + prediction: Predictions, + *, + prefix_last: object, + boundary: object, + fold: Mapping[str, Any], + specimen: Mapping[str, Any], + protocol: Mapping[str, Any], +) -> dict[str, tuple[np.ndarray, np.ndarray]]: + """Construct every registered arm without reading a target outcome.""" + validate_protocol(protocol) + residuals = _mapping(fold["source_residual_variance_m2"], name="residuals") + nominal_mean = _finite_array(prediction.nominal, name="prediction.nominal", ndim=3) + if nominal_mean.shape != prediction.bank.shape[1:]: + raise ValueError("nominal trajectory and model bank disagree") + result: dict[str, tuple[np.ndarray, np.ndarray]] = { + "nominal_physics": ( + nominal_mean, + _variance_field( + prediction.inputs, + nominal_mean.shape, + residuals["nominal_physics"], + ), + ) + } + last_mean = _last_residual_mean(prediction, prefix_last, boundary) + result["last_residual"] = ( + last_mean, + _variance_field(prediction.inputs, last_mean.shape, residuals["last_residual"]), + ) + bayesian_residual = residuals["bayesian"] + single = _mapping(specimen["single_probe_weights"], name="single_probe_weights") + for condition in protocol["probe_conditions"]: + if condition not in single: + raise ValueError(f"missing single-probe weights for {condition}") + result[f"single_probe_{condition}"] = weighted_belief( + prediction, single[condition], bayesian_residual + ) + states = _mapping(specimen["policy_states"], name="policy_states") + for policy in protocol["probe_policies"]: + policy_states = _mapping(states[policy], name=f"policy_states[{policy}]") + for budget in protocol["probe_budgets"]: + state = _mapping( + policy_states[str(budget)], + name=f"policy_states[{policy}][{budget}]", + ) + selected = tuple(str(value) for value in state["selected_actions"]) + if len(selected) != int(budget) or len(set(selected)) != len(selected): + raise ValueError("selected action roster does not match its budget") + if not set(selected).issubset(protocol["probe_conditions"]): + raise ValueError("policy state contains an unregistered action") + result[f"{policy}_k{budget}"] = weighted_belief( + prediction, state["weights"], bayesian_residual + ) + expected = set(arm_specs(protocol)) + if set(result) != expected: + raise ValueError("constructed belief roster differs from the protocol") + for mean, variance in result.values(): + if mean.shape != nominal_mean.shape or variance.shape != nominal_mean.shape: + raise ValueError("belief trajectory shape mismatch") + if not np.all(np.isfinite(mean)) or not np.all(np.isfinite(variance)): + raise ValueError("belief trajectory contains nonfinite values") + if np.any(variance <= 0.0): + raise ValueError("belief variance must be positive") + return result + + +__all__ = [ + "active_mask", + "arm_specs", + "array_digest", + "belief_digest", + "build_belief_arms", + "calibrated_residuals", + "loss_vector", + "posterior_temperature", + "validate_protocol", + "weighted_belief", +] diff --git a/experiments/tracking_cloth_deformation_v1/data.py b/experiments/tracking_cloth_deformation_v1/data.py new file mode 100644 index 000000000..0ac4ab4de --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/data.py @@ -0,0 +1,368 @@ +"""Read-only archive audit and causal CSV views for Zenodo 14644526. + +No MATLAB code is executed. The parser follows the paper's Frame, Time, XYZ +layout. A target input view converts future corner columns only; future free +marker values are converted exclusively by the post-seal scoring reader. +""" + +from __future__ import annotations + +import csv +import hashlib +import itertools +import json +import re +import zipfile +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +import numpy as np + +FREE_NAME = re.compile( + r"^(cotton|denim|polyester|wool)_(A2|A3)_(shake|twist)_" + r"(fast|slow)_(hands|hanger)\.csv$", + re.IGNORECASE, +) + + +@dataclass(frozen=True) +class Case: + path: Path + material: str + size: str + motion: str + speed: str + grasp: str + + @property + def specimen(self) -> str: + return f"{self.material}_{self.size}" + + @property + def condition(self) -> str: + return f"{self.speed}_{self.grasp}" + + @property + def markers(self) -> int: + return 20 if self.size == "A2" else 12 + + +def digest(path: Path, algorithm: str = "sha256") -> str: + # MD5 is required only to match the publisher's archive checksum. + h = hashlib.new(algorithm, usedforsecurity=False) + with path.open("rb") as stream: + for block in iter(lambda: stream.read(1024 * 1024), b""): + h.update(block) + return h.hexdigest() + + +def object_digest(value: Any) -> str: + return hashlib.sha256( + json.dumps( + value, sort_keys=True, separators=(",", ":"), allow_nan=False + ).encode() + ).hexdigest() + + +def write_json(path: Path, value: Any) -> None: + path.write_text(json.dumps(value, indent=2, sort_keys=True, allow_nan=False) + "\n") + + +def audit_dataset( + root: Path, protocol: dict[str, Any] +) -> tuple[list[Case], dict[str, Any]]: + """Hash bytes, verify extraction against ZIP, and freeze the filename roster. + + Hashing/ZIP integrity reads bytes but does not interpret target measurements. + Files are never changed, downloaded, extracted, or chmod'ed by this runner. + """ + root = root.resolve(strict=True) + files = sorted(p for p in root.rglob("*") if p.is_file()) + for path in files: + if not path.resolve().is_relative_to(root): + raise ValueError(f"Dataset entry escapes the cache: {path.name}") + csvs = [p for p in files if p.suffix.lower() == ".csv"] + if len(csvs) != protocol["csv_count"]: + raise ValueError(f"Expected {protocol['csv_count']} CSVs, found {len(csvs)}") + if len({p.name.lower() for p in csvs}) != len(csvs): + raise ValueError("Ambiguous duplicate CSV basenames") + archives = [p for p in files if p.suffix.lower() == ".zip"] + matching = [p for p in archives if digest(p, "md5") == protocol["archive_md5"]] + if len(matching) != 1: + raise ValueError("Expected one retained ZIP matching the published MD5") + archive = matching[0] + readers = [p for p in files if p.name.lower() == "read_data.m"] + licenses = [p for p in files if p.name.lower() == "license.txt"] + if len(readers) != 4 or len(licenses) != 1: + raise ValueError("Expected four read_data.m readers and one License.txt") + license_text = licenses[0].read_text(encoding="utf-8-sig") + normalized = license_text.lower().replace(" ", "") + if not any( + s in normalized for s in ("by-nc-sa", "noncommercial", "non-commercial") + ): + raise ValueError( + "Included license differs from the declared noncommercial policy" + ) + hashes = {p.name.lower(): digest(p) for p in csvs} + with zipfile.ZipFile(archive) as zipped: + bad = zipped.testzip() + if bad: + raise ValueError(f"ZIP integrity failure: {bad}") + entries = [z for z in zipped.infolist() if z.filename.lower().endswith(".csv")] + if len(entries) != len(csvs): + raise ValueError("ZIP and extracted CSV inventories disagree") + seen = set() + for entry in entries: + name = Path(entry.filename).name.lower() + if name in seen or name not in hashes: + raise ValueError("Ambiguous ZIP CSV identity") + seen.add(name) + if hashlib.sha256(zipped.read(entry)).hexdigest() != hashes[name]: + raise ValueError(f"Extracted bytes differ from verified ZIP: {name}") + cases = [] + for path in csvs: + match = FREE_NAME.fullmatch(path.name) + if match: + material, size, motion, speed, grasp = match.groups() + cases.append( + Case( + path, + material.lower(), + size.upper(), + motion.lower(), + speed.lower(), + grasp.lower(), + ) + ) + expected = set( + itertools.product( + protocol["materials"], + protocol["sizes"], + ["shake", "twist"], + protocol["speeds"], + protocol["grasps"], + ) + ) + actual = {(c.material, c.size, c.motion, c.speed, c.grasp) for c in cases} + if actual != expected or len(cases) != 64: + raise ValueError("The complete 64-recording free-hanging factorial is required") + inventory = { + "dataset_record": protocol["dataset_record"], + "archive_name": archive.name, + "archive_md5": protocol["archive_md5"], + "archive_sha256": digest(archive), + "csv_count": len(csvs), + "source_count": 32, + "target_count": 32, + "unused_count": len(csvs) - 64, + "csv_sha256": hashes, + "reader_sha256": {str(p.relative_to(root)): digest(p) for p in readers}, + "license_sha256": digest(licenses[0]), + "included_license_text": license_text, + "license_policy": protocol["license_policy"], + "target_numeric_outcomes_read": False, + } + inventory["inventory_id"] = object_digest(inventory) + return sorted(cases, key=lambda c: c.path.name), inventory + + +def _rows(path: Path, markers: int): + """Yield the documented numeric Frame, Time, XYZ rows, preserving blanks.""" + width = 2 + 3 * markers + started = False + last_t = -np.inf + last_frame = -np.inf + with path.open(encoding="utf-8-sig", newline="") as stream: + for row in csv.reader(stream): + if not row or not any(cell.strip() for cell in row): + continue + try: + frame, time = float(row[0]), float(row[1]) + except (ValueError, IndexError): + if started: + raise ValueError( + f"Nonnumeric row after data start: {path.name}" + ) from None + continue + started = True + if ( + not np.isfinite([frame, time]).all() + or frame != int(frame) + or frame <= last_frame + or time <= last_t + ): + raise ValueError(f"Invalid frame/time order: {path.name}") + if len(row) < width or any(cell.strip() for cell in row[width:]): + raise ValueError( + f"Expected Frame, Time and {markers} XYZ triplets: {path.name}" + ) + last_t, last_frame = time, frame + yield time, row[2:width] + if not started: + raise ValueError(f"No numeric rows in {path.name}") + + +def _positions(cells: list[str], indices: np.ndarray) -> np.ndarray: + values = [] + for marker in indices: + triple = [ + float(cells[3 * int(marker) + d]) + if cells[3 * int(marker) + d].strip() + else np.nan + for d in range(3) + ] + if np.isinf(triple).any(): + raise ValueError("Infinite coordinate") + values.append(triple if np.isfinite(triple).all() else [np.nan] * 3) + return np.asarray(values, dtype=float) + + +def read_prefix(case: Case, seconds: float) -> tuple[np.ndarray, np.ndarray]: + times, values = [], [] + start = None + for time, cells in _rows(case.path, case.markers): + if start is None: + start = time + if time - start > seconds + 1e-8: + break + times.append(time) + values.append(_positions(cells, np.arange(case.markers))) + return np.asarray(times), np.asarray(values) + + +def infer_source_scale(case: Case, positions: np.ndarray) -> float: + """Resolve m/cm/mm from source-only initial geometry; never from a target.""" + complete = np.flatnonzero(np.isfinite(positions).all(axis=(1, 2))) + if not len(complete): + raise ValueError(f"No complete initialization frame: {case.path.name}") + first = positions[complete[0]] + diameter = np.linalg.norm(first[:, None] - first[None, :], axis=2).max() + nominal = np.hypot(0.42, 0.594) if case.size == "A2" else np.hypot(0.297, 0.42) + allowed = [s for s in (1.0, 0.01, 0.001) if 0.4 < diameter * s / nominal < 1.6] + if len(allowed) != 1: + raise ValueError( + "Ambiguous coordinate units; inspect source readers before revising protocol" + ) + return allowed[0] + + +def layout(initial: np.ndarray, size: str) -> tuple[np.ndarray, np.ndarray]: + """Initial-frame grid only: four/five descending-z rows, three/four columns. + + A deliberately explicit assumption of the first pilot, not claimed to be the + official MATLAB marker-ordering implementation. Ambiguous geometry fails. + """ + rows, cols = (5, 4) if size == "A2" else (4, 3) + xy = initial[:, :2] - initial[:, :2].mean(axis=0) + _, _, vh = np.linalg.svd(xy, full_matrices=False) + horizontal = vh[0] + horizontal *= 1 if horizontal[np.argmax(abs(horizontal))] >= 0 else -1 + levels = np.argsort(-initial[:, 2], kind="stable").reshape(rows, cols) + order = np.concatenate( + [row[np.argsort(xy[row] @ horizontal, kind="stable")] for row in levels] + ) + grid = initial[order].reshape(rows, cols, 3) + vertical_steps = grid[:-1, :, 2].mean(axis=1) - grid[1:, :, 2].mean(axis=1) + horizontal_steps = np.linalg.norm(np.diff(grid, axis=1), axis=2) + if ( + np.min(vertical_steps) < 0.025 + or np.min(horizontal_steps) < 0.025 + or np.max(horizontal_steps) > 0.25 + ): + raise ValueError("Initial grid/corner assignment is unsupported") + # This pilot assumes a near-vertical regular initial mesh, not arbitrary folds. + if np.max(np.ptp(grid[:, :, 2], axis=1)) > 0.75 * np.median(vertical_steps): + raise ValueError( + "Initial row ordering is ambiguous; no outcome-based remapping" + ) + return order, np.array([0, cols - 1], dtype=int) + + +@dataclass(frozen=True) +class Inputs: + times: np.ndarray + prefix: np.ndarray + boundary: np.ndarray + order: np.ndarray + corners: np.ndarray + cutoff: int + initial_time: float + scale: float + + +def input_view(case: Case, protocol: dict[str, Any], scale: float) -> Inputs: + """Only prefix XYZ and future prescribed-corner XYZ enter this view.""" + seconds = protocol["prefix_seconds"] + time0, pos0 = read_prefix(case, seconds) + complete = np.flatnonzero(np.isfinite(pos0).all(axis=(1, 2))) + if not len(complete): + raise ValueError(f"No complete prefix frame: {case.path.name}") + first = int(complete[0]) + if time0[first] - time0[0] > protocol["initial_complete_frame_deadline_seconds"]: + raise ValueError("Late complete initialization frame") + initial_time = float(time0[first]) + order, corners = layout(pos0[first] * scale, case.size) + start_time = float(time0[0]) + stride = protocol["sample_stride"] + end = start_time + seconds + protocol["forecast_seconds"] + times, prefix, boundary = [], [], [] + last_boundary = None + row_index = 0 + for time, cells in _rows(case.path, case.markers): + if time < initial_time - 1e-9: + continue + if time > end + 1e-8: + break + use = row_index % stride == 0 + row_index += 1 + if not use: + continue + b = _positions(cells, order[corners]) * scale + if last_boundary is not None: + b = np.where(np.isfinite(b), b, last_boundary) + if not np.isfinite(b).all(): + raise ValueError("Missing initial driven corner") + last_boundary = b.copy() + times.append(time) + boundary.append(b) + if time <= start_time + seconds + 1e-8: + p = _positions(cells, order) * scale + if prefix: + p = np.where(np.isfinite(p), p, prefix[-1]) + if not np.isfinite(p).all(): + raise ValueError("Nonfinite causal initialization") + prefix.append(p) + times_array = np.asarray(times) + if len(prefix) < 5 or len(times) <= len(prefix) + 5: + raise ValueError("Insufficient prefix or forecast frames") + dt = np.diff(times_array) + if not np.allclose(dt, stride / 120.0, rtol=0.05, atol=1e-4): + raise ValueError("Sampling does not match the frozen 120 Hz/stride contract") + if times_array[-1] < end - 2 * stride / 120.0: + raise ValueError("Recording does not cover the complete frozen horizon") + return Inputs( + times_array, + np.asarray(prefix), + np.asarray(boundary), + order, + corners, + len(prefix) - 1, + initial_time, + scale, + ) + + +def scoring_view(case: Case, inputs: Inputs) -> np.ndarray: + """Called only by source training or after the complete target prediction seal.""" + rows = [] + index = 0 + for time, cells in _rows(case.path, case.markers): + if index == len(inputs.times): + break + if abs(time - inputs.times[index]) <= 1e-7: + rows.append(_positions(cells, inputs.order) * inputs.scale) + index += 1 + if index != len(inputs.times): + raise ValueError("Scoring timestamps do not reproduce the sealed input view") + return np.asarray(rows) diff --git a/experiments/tracking_cloth_deformation_v1/model.py b/experiments/tracking_cloth_deformation_v1/model.py new file mode 100644 index 000000000..21b33e36d --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/model.py @@ -0,0 +1,301 @@ +"""Small NumPy spring-mesh pilot and source-only generalized Bayes. + +This is an explicitly limited equal-marker-mass spring model, not PhysTwin, +clothilde-sim, a FEM reproduction, or a newly validated material estimator. +""" + +from __future__ import annotations + +import itertools +from dataclasses import dataclass +from typing import Any + +import numpy as np + +from .data import Inputs + +ARMS = ( + "persistence", + "nominal_physics", + "last_residual", + "nominal_state_injection", + "map_physics", + "bayesian_physics", + "guarded_bayesian_physics", +) + + +def parameter_bank(protocol: dict[str, Any]) -> list[tuple[float, float]]: + return list( + itertools.product(protocol["stiffness_per_mass"], protocol["damping_per_mass"]) + ) + + +def edges(markers: int) -> tuple[np.ndarray, np.ndarray]: + rows, cols = (5, 4) if markers == 20 else (4, 3) + links, weights = [], [] + for r in range(rows): + for c in range(cols): + for dr, dc, weight in ( + (0, 1, 1.0), + (1, 0, 1.0), + (1, 1, 0.5), + (1, -1, 0.5), + (0, 2, 0.1), + (2, 0, 0.1), + ): + if 0 <= r + dr < rows and 0 <= c + dc < cols: + links.append((r * cols + c, (r + dr) * cols + c + dc)) + weights.append(weight) + return np.asarray(links, dtype=int), np.asarray(weights) + + +def velocity(times: np.ndarray, positions: np.ndarray) -> np.ndarray: + t = times - times.mean() + return np.einsum("t,tnd->nd", t, positions) / np.dot(t, t) + + +def rollout( + inputs: Inputs, + parameters: tuple[float, float], + protocol: dict[str, Any], + inject: bool, +) -> np.ndarray: + """Symplectic spring rollout with recorded Dirichlet corner input. + + Rest edge lengths use the initial frame, not the forecast outcomes. A state + injection uses the final permitted prefix position and backward-only velocity. + """ + links, relative_k = edges(len(inputs.order)) + left, right = links.T + rest = np.linalg.norm(inputs.prefix[0, right] - inputs.prefix[0, left], axis=1) + if np.min(rest) <= 1e-6: + raise ValueError("Degenerate initial spring") + k, damping = parameters + start = inputs.cutoff if inject else 0 + x = inputs.prefix[start].copy() + if inject: + v = velocity( + inputs.times[start - 4 : start + 1], inputs.prefix[start - 4 : start + 1] + ) + else: + v = velocity(inputs.times[:5], inputs.prefix[:5]) + result = np.empty((len(inputs.times), len(inputs.order), 3)) + result[: start + 1] = inputs.prefix[: start + 1] + origin = inputs.prefix[0].mean(axis=0) + substeps = protocol["integration_substeps"] + for t in range(start + 1, len(inputs.times)): + full_dt = inputs.times[t] - inputs.times[t - 1] + dt = full_dt / substeps + boundary_v = (inputs.boundary[t] - inputs.boundary[t - 1]) / full_dt + for sub in range(1, substeps + 1): + delta = x[right] - x[left] + lengths = np.linalg.norm(delta, axis=1) + force = (k * relative_k * (lengths - rest) / np.maximum(lengths, 1e-9))[ + :, None + ] * delta + acceleration = -damping * v + acceleration[:, 2] -= protocol["gravity_m_s2"] + np.add.at(acceleration, left, force) + np.add.at(acceleration, right, -force) + v += dt * acceleration + x += dt * v + fraction = sub / substeps + x[inputs.corners] = (1 - fraction) * inputs.boundary[ + t - 1 + ] + fraction * inputs.boundary[t] + v[inputs.corners] = boundary_v + if not np.isfinite(x).all() or np.max(np.linalg.norm(x - origin, axis=1)) > 10: + raise ValueError("Numerically invalid rollout; no silent case deletion") + result[t] = x + return result + + +@dataclass(frozen=True) +class Predictions: + inputs: Inputs + nominal: np.ndarray + bank: np.ndarray + + +def predict(inputs: Inputs, protocol: dict[str, Any]) -> Predictions: + nominal = rollout(inputs, tuple(protocol["nominal_parameters"]), protocol, False) + bank = np.stack( + [ + rollout(inputs, parameters, protocol, True) + for parameters in parameter_bank(protocol) + ] + ) + return Predictions(inputs, nominal, bank) + + +def masks(inputs: Inputs, truth: np.ndarray) -> np.ndarray: + valid = np.isfinite(truth).all(axis=2) + valid[: inputs.cutoff + 1] = False + valid[:, inputs.corners] = False + if not np.any(valid): + raise ValueError("No evaluable free-marker forecast samples") + return valid + + +def squared_error(mean: np.ndarray, truth: np.ndarray, valid: np.ndarray) -> float: + if not np.isfinite(mean).all(): + raise ValueError("Nonfinite prediction") + return float(np.mean(np.sum((mean[valid] - truth[valid]) ** 2, axis=1))) + + +def source_weights(losses: np.ndarray, floor_m: float) -> np.ndarray: + """Gibbs update: one normalized trajectory loss per source recording. + + The source-derived temperature is fixed before target use. Marker/time rows + are not counted as independent likelihood groups or inferential replicates. + """ + temperature = max(float(np.median(np.min(losses, axis=1))), floor_m**2) + logits = -np.sum(losses, axis=0) / (2 * temperature) + weights = np.exp(logits - np.max(logits)) + return weights / weights.sum() + + +def means( + predictions: Predictions, weights: np.ndarray, protocol: dict[str, Any] +) -> dict[str, np.ndarray]: + inputs = predictions.inputs + residual = inputs.prefix[-1] - predictions.nominal[inputs.cutoff] + nominal_index = parameter_bank(protocol).index( + tuple(protocol["nominal_parameters"]) + ) + result = { + "persistence": np.broadcast_to( + inputs.prefix[-1], predictions.nominal.shape + ).copy(), + "nominal_physics": predictions.nominal, + "last_residual": predictions.nominal + residual, + "nominal_state_injection": predictions.bank[nominal_index], + "map_physics": predictions.bank[int(np.argmax(weights))], + "bayesian_physics": np.einsum("k,ktnd->tnd", weights, predictions.bank), + } + # The boundary is a conditioning input, never a scored output. + for name in ("persistence", "last_residual"): + result[name][:, inputs.corners] = inputs.boundary + return result + + +def horizon_bins(inputs: Inputs) -> np.ndarray: + duration = inputs.times[-1] - inputs.times[inputs.cutoff] + return np.clip( + ((inputs.times - inputs.times[inputs.cutoff]) / duration * 3).astype(int), 0, 2 + ) + + +def fit_specimen( + records: list[tuple[Predictions, np.ndarray]], protocol: dict[str, Any] +) -> dict[str, Any]: + """Leave one speed/grasp source recording out; never access target outcomes.""" + if len(records) != 4: + raise ValueError( + "Exactly four shaking source conditions are required per specimen" + ) + loss = np.asarray( + [ + [squared_error(p, truth, masks(pred.inputs, truth)) for p in pred.bank] + for pred, truth in records + ] + ) + oof = [] + residual_squares = {arm: [[], [], []] for arm in ARMS[:-1]} + for held, (prediction, truth) in enumerate(records): + w = source_weights( + np.delete(loss, held, axis=0), protocol["measurement_floor_m"] + ) + arm_means = means(prediction, w, protocol) + valid = masks(prediction.inputs, truth) + bins = horizon_bins(prediction.inputs) + oof.append( + { + arm: np.sqrt(squared_error(mean, truth, valid)) + for arm, mean in arm_means.items() + } + ) + ensemble_var = np.einsum( + "k,ktnd->tnd", w, (prediction.bank - arm_means["bayesian_physics"]) ** 2 + ) + for arm, mean in arm_means.items(): + for b in range(3): + select = valid & (bins[:, None] == b) + if not np.any(select): + raise ValueError("Empty source calibration horizon bin") + error2 = (mean[select] - truth[select]) ** 2 + if arm == "bayesian_physics": + error2 = error2 - ensemble_var[select] + # Equal recording contribution to each calibration bin. + residual_squares[arm][b].append(float(np.mean(error2))) + baseline = np.asarray([row["nominal_physics"] for row in oof]) + residual = np.asarray([row["last_residual"] for row in oof]) + candidate = np.asarray([row["bayesian_physics"] for row in oof]) + reference = min(float(baseline.mean()), float(residual.mean())) + accepted = bool( + np.all(candidate <= baseline) + and candidate.mean() < (1 - protocol["guard_minimum_relative_gain"]) * reference + ) + noise = { + arm: [ + max(float(np.mean(values)), protocol["measurement_floor_m"] ** 2) + for values in bins + ] + for arm, bins in residual_squares.items() + } + return { + "source_posterior_weights": source_weights( + loss, protocol["measurement_floor_m"] + ).tolist(), + "guard_accepts": accepted, + "guard_basis": "OOF mean beats nominal and last_residual by frozen margin; no OOF nominal regression", + "oof_record_rmse_m": oof, + "source_residual_variance_m2": noise, + } + + +def complete_beliefs( + prediction: Predictions, fit: dict[str, Any], protocol: dict[str, Any] +) -> dict[str, tuple[np.ndarray, np.ndarray]]: + weights = np.asarray(fit["source_posterior_weights"]) + arm_means = means(prediction, weights, protocol) + bins = horizon_bins(prediction.inputs) + beliefs = {} + for arm, mean in arm_means.items(): + variance = np.broadcast_to( + np.asarray(fit["source_residual_variance_m2"][arm])[bins, None, None], + mean.shape, + ).copy() + if arm == "bayesian_physics": + variance += np.einsum("k,ktnd->tnd", weights, (prediction.bank - mean) ** 2) + beliefs[arm] = (mean, variance) + # Preserve both mean and covariance together; do not reconstruct on rejection. + chosen = "bayesian_physics" if fit["guard_accepts"] else "nominal_physics" + beliefs["guarded_bayesian_physics"] = beliefs[chosen] + return beliefs + + +def score( + mean: np.ndarray, variance: np.ndarray, truth: np.ndarray, inputs: Inputs +) -> dict[str, float | int]: + valid = masks(inputs, truth) + e = mean[valid] - truth[valid] + var = variance[valid] + if not np.isfinite(var).all() or np.any(var <= 0): + raise ValueError("Invalid predictive variance") + return { + "rmse_mm": 1000 * float(np.sqrt(np.mean(np.sum(e**2, axis=1)))), + "mean_marker_error_mm": 1000 * float(np.mean(np.linalg.norm(e, axis=1))), + "coordinate_nll": float(np.mean(0.5 * (np.log(2 * np.pi * var) + e**2 / var))), + "coordinate_90_coverage": float( + np.mean(abs(e) <= 1.6448536269514722 * np.sqrt(var)) + ), + "mean_full_90_width_mm": 1000 + * float(np.mean(2 * 1.6448536269514722 * np.sqrt(var))), + "free_marker_samples": int(valid.sum()), + "missing_free_marker_samples": int( + (len(inputs.times) - inputs.cutoff - 1) * (len(inputs.order) - 2) + - valid.sum() + ), + } diff --git a/experiments/tracking_cloth_deformation_v1/protocol.json b/experiments/tracking_cloth_deformation_v1/protocol.json new file mode 100644 index 000000000..49dd4a15d --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/protocol.json @@ -0,0 +1,30 @@ +{ + "study_id": "tracking-cloth-shake-to-twist-pilot-v1", + "dataset_record": "14644526", + "archive_md5": "b4868b702f8a42b2ea1069d0f1a3b8f6", + "csv_count": 120, + "source_motion": "shake", + "target_motion": "twist", + "materials": ["cotton", "denim", "polyester", "wool"], + "sizes": ["A2", "A3"], + "speeds": ["fast", "slow"], + "grasps": ["hands", "hanger"], + "sample_stride": 4, + "prefix_seconds": 1.0, + "forecast_seconds": 5.0, + "initial_complete_frame_deadline_seconds": 0.25, + "stiffness_per_mass": [100.0, 400.0, 1600.0], + "damping_per_mass": [0.5, 2.0, 8.0], + "nominal_parameters": [400.0, 2.0], + "integration_substeps": 8, + "gravity_m_s2": 9.80665, + "guard_minimum_relative_gain": 0.01, + "measurement_floor_m": 0.001, + "bootstrap_repetitions": 10000, + "bootstrap_seed": 14644526, + "information_regime": "measured-corner-boundary-conditioned; all-marker initialization prefix only", + "evidence_class": "public-real-data-pilot; prior exposure not assumed absent", + "license_policy": "CC-BY-NC-SA-4.0 pending author clarification of conflicting CC-BY-SA-4.0 metadata", + "raw_data_upload": false, + "paper_claim_authorized": false +} diff --git a/experiments/tracking_cloth_deformation_v1/requirements.txt b/experiments/tracking_cloth_deformation_v1/requirements.txt new file mode 100644 index 000000000..fbe88a68b --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/requirements.txt @@ -0,0 +1,2 @@ +# Isolated CPU pilot; no CUDA, MATLAB, or system-wide installation required. +numpy==1.26.4 diff --git a/experiments/tracking_cloth_deformation_v1/run.py b/experiments/tracking_cloth_deformation_v1/run.py new file mode 100644 index 000000000..e0d5b088c --- /dev/null +++ b/experiments/tracking_cloth_deformation_v1/run.py @@ -0,0 +1,503 @@ +"""Maintained CLI for read-only inventory, source fit, prediction seal and scoring.""" + +from __future__ import annotations + +import argparse +import csv +import json +import os +import platform +import sys +import traceback +from concurrent.futures import ProcessPoolExecutor +from datetime import datetime, timezone +from pathlib import Path +from typing import Any + +import numpy as np + +from .data import ( + Inputs, + audit_dataset, + digest, + infer_source_scale, + input_view, + object_digest, + read_prefix, + scoring_view, + write_json, +) +from .model import ARMS, complete_beliefs, fit_specimen, predict, score + +HERE = Path(__file__).resolve().parent +METRICS = ( + "rmse_mm", + "mean_marker_error_mm", + "coordinate_nll", + "coordinate_90_coverage", + "mean_full_90_width_mm", +) + + +def now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def implementation() -> dict[str, str]: + return {p.name: digest(p) for p in sorted(HERE.glob("*.py"))} + + +def source_record(args): + case, protocol, scale = args + inputs = input_view(case, protocol, scale) + prediction = predict(inputs, protocol) + return case.specimen, case.path.name, prediction, scoring_view(case, inputs) + + +def save_csv(path: Path, rows: list[dict[str, Any]]) -> None: + if not rows: + raise ValueError("Refusing an empty result table") + with path.open("w", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=list(rows[0])) + writer.writeheader() + writer.writerows(rows) + + +def prepare( + root: Path, output: Path, protocol: dict[str, Any], stage: str, workers: int = 1 +) -> None: + if stage not in ("inventory", "source", "predict"): + raise ValueError("Unknown preparation stage") + root = root.resolve(strict=True) + output = output.resolve() + if output.is_relative_to(root) or root.is_relative_to(output): + raise ValueError("Output and dataset must be disjoint directory trees") + output.mkdir(parents=True, exist_ok=False) + write_json(output / "protocol.json", protocol) + provenance = { + "created_at": now(), + "protocol_id": object_digest(protocol), + "implementation_sha256": implementation(), + "python": sys.version, + "numpy": np.__version__, + "platform": platform.platform(), + "github_sha": os.environ.get("GITHUB_SHA"), + "github_run_id": os.environ.get("GITHUB_RUN_ID"), + "github_run_attempt": os.environ.get("GITHUB_RUN_ATTEMPT"), + "runner_name": os.environ.get("RUNNER_NAME"), + "target_numeric_outcomes_read": False, + "evidence_class": protocol["evidence_class"], + "paper_claim_authorized": False, + } + write_json(output / "run_manifest.json", provenance) + cases, inventory = audit_dataset(root, protocol) + write_json(output / "dataset_manifest.json", inventory) + (output / "DATA_LICENSE.txt").write_text(inventory["included_license_text"]) + report = [ + "# Tracking Cloth Deformation: public-data pilot", + "", + f"Study: `{protocol['study_id']}`", + "", + "120 verified CSVs; 32 shaking source recordings, 32 twisting targets,", + "and 56 collision recordings reserved and not numerically read.", + "", + "Dataset cache is read-only. Archive/extracted-byte hashing is not", + "numeric outcome evaluation. Included noncommercial license governs", + "this run pending author clarification; raw recordings are not uploaded.", + "", + "This is a reduced spring-mesh pilot, not a PhysTwin/FEM reproduction.", + "No new acquisition or paper claim is created.", + "", + ] + (output / "report.md").write_text("\n".join(report)) + if stage == "inventory": + return + source = [c for c in cases if c.motion == "shake"] + scales = [ + infer_source_scale(c, read_prefix(c, protocol["prefix_seconds"])[1]) + for c in source + ] + if len(set(scales)) != 1: + raise ValueError("Source recordings disagree about metric coordinate units") + scale = scales[0] + tasks = [(c, protocol, scale) for c in source] + if workers > 1: + with ProcessPoolExecutor(max_workers=workers) as pool: + records = list(pool.map(source_record, tasks)) + else: + records = [source_record(task) for task in tasks] + fitted = {} + source_rows = [] + for specimen in sorted({c.specimen for c in source}): + subset = [ + (name, pred, truth) + for group, name, pred, truth in records + if group == specimen + ] + fitted[specimen] = fit_specimen( + [(pred, truth) for _, pred, truth in subset], protocol + ) + fitted[specimen]["source_recordings"] = [name for name, _, _ in subset] + for name, oof in zip( + fitted[specimen]["source_recordings"], + fitted[specimen]["oof_record_rmse_m"], + strict=True, + ): + source_rows.extend( + { + "recording": name, + "specimen": specimen, + "arm": arm, + "oof_rmse_mm": 1000 * float(value), + } + for arm, value in oof.items() + ) + freeze = { + "protocol_id": object_digest(protocol), + "inventory_id": inventory["inventory_id"], + "implementation_sha256": implementation(), + "coordinate_scale_to_m": scale, + "fitted_at": now(), + "specimens": fitted, + "target_outcomes_used": False, + "guard_is_empirical_source_rule_not_safety_certificate": True, + } + write_json(output / "source_fit.json", freeze) + save_csv(output / "source_scores.csv", source_rows) + accepted = sum(int(f["guard_accepts"]) for f in fitted.values()) + report.extend( + [ + "## Source-only qualification", + "", + f"Coordinate scale to metres: `{scale}` (inferred from source initialization only).", + f"Empirical source guard accepts {accepted}/8 specimen candidates.", + "Each specimen uses four-fold leave-one-speed/grasp-recording-out fitting.", + "These folds select the model/guard; they are not independent confirmation.", + "No twisting free-marker forecast outcome has been evaluated.", + "", + ] + ) + (output / "report.md").write_text("\n".join(report)) + if stage == "source": + return + private = output / "private_predictions" + private.mkdir(mode=0o700) + predictions = {} + for case in (c for c in cases if c.motion == "twist"): + inputs = input_view(case, protocol, scale) + beliefs = complete_beliefs( + predict(inputs, protocol), fitted[case.specimen], protocol + ) + arrays = {f"{arm}_mean": beliefs[arm][0] for arm in ARMS} + arrays.update({f"{arm}_variance": beliefs[arm][1] for arm in ARMS}) + arrays.update( + { + "times": inputs.times, + "order": inputs.order, + "corners": inputs.corners, + "cutoff": np.array(inputs.cutoff), + "scale": np.array(scale), + } + ) + artifact = private / f"{case.path.stem}.npz" + np.savez_compressed(artifact, **arrays) + predictions[case.path.name] = { + "artifact": str(artifact.relative_to(output)), + "sha256": digest(artifact), + "specimen": case.specimen, + "guard_accepts": fitted[case.specimen]["guard_accepts"], + "corner_raw_column_indices": inputs.order[inputs.corners].tolist(), + "causal_cutoff_seconds": float(inputs.times[inputs.cutoff]), + } + if len(predictions) != 32: + raise ValueError("Refusing an incomplete target prediction seal") + seal = { + "sealed_at": now(), + "protocol_id": object_digest(protocol), + "inventory_id": inventory["inventory_id"], + "source_fit_sha256": digest(output / "source_fit.json"), + "implementation_sha256": implementation(), + "predictions": predictions, + "future_free_marker_outcomes_read": False, + "future_driven_corner_coordinates_used": True, + "initialization_prefix_all_markers_used": True, + "prior_public_outcome_exposure": "unknown; no fresh-confirmation claim", + } + write_json(output / "prediction_seal.json", seal) + report.extend( + [ + "## Predictions sealed", + "", + "All 32 target batches are sealed before scoring.", + "Only timestamps, the initialization prefix and future prescribed corners", + "entered prediction. Forecasts stay local and are not in the upload bundle.", + "", + ] + ) + (output / "report.md").write_text("\n".join(report)) + + +def aggregate( + rows: list[dict[str, Any]], protocol: dict[str, Any] +) -> tuple[list[dict[str, Any]], dict[str, Any]]: + specimens = sorted({row["specimen"] for row in rows}) + if len(specimens) != 8 or len(rows) != 32 * len(ARMS): + raise ValueError("Incomplete roster; no pooled partial result is authorized") + table = [] + for specimen in specimens: + for arm in ARMS: + subset = [ + row for row in rows if row["specimen"] == specimen and row["arm"] == arm + ] + if len(subset) != 4 or len({row["recording"] for row in subset}) != 4: + raise ValueError("Missing or duplicate speed/grasp condition") + table.append( + { + "specimen": specimen, + "arm": arm, + **{ + metric: float(np.mean([row[metric] for row in subset])) + for metric in METRICS + }, + } + ) + summary = { + arm: { + metric: float(np.mean([row[metric] for row in table if row["arm"] == arm])) + for metric in METRICS + } + for arm in ARMS + } + rng = np.random.default_rng(protocol["bootstrap_seed"]) + resamples = rng.integers(0, 8, size=(protocol["bootstrap_repetitions"], 8)) + contrasts = {} + for comparator in ("nominal_physics", "last_residual", "map_physics"): + diffs = np.array( + [ + next( + r["rmse_mm"] + for r in table + if r["specimen"] == s and r["arm"] == "guarded_bayesian_physics" + ) + - next( + r["rmse_mm"] + for r in table + if r["specimen"] == s and r["arm"] == comparator + ) + for s in specimens + ] + ) + material_diffs = np.array( + [ + np.mean( + [ + diffs[i] + for i, s in enumerate(specimens) + if s.startswith(material + "_") + ] + ) + for material in protocol["materials"] + ] + ) + material_samples = rng.integers( + 0, 4, size=(protocol["bootstrap_repetitions"], 4) + ) + contrasts[comparator] = { + "guarded_minus_comparator_rmse_mm": float(diffs.mean()), + "specimen_bootstrap_95_interval_mm": np.quantile( + diffs[resamples].mean(axis=1), [0.025, 0.975] + ).tolist(), + "material_cluster_sensitivity_95_interval_mm": np.quantile( + material_diffs[material_samples].mean(axis=1), [0.025, 0.975] + ).tolist(), + "specimen_wins": int((diffs < 0).sum()), + "specimen_ties": int((diffs == 0).sum()), + "specimen_losses": int((diffs > 0).sum()), + "worst_specimen_regret_mm": float(diffs.max()), + } + return table, { + "arms": summary, + "contrasts": contrasts, + "inferential_unit": "8 material-size specimens; 4-material sensitivity also reported", + "interval_interpretation": "exploratory paired percentile bootstrap; not simultaneous; small cluster counts", + "aggregation": "equal recordings within specimen, then equal specimens; no frame pseudoreplication", + } + + +def score_run(root: Path, output: Path) -> None: + root, output = root.resolve(strict=True), output.resolve(strict=True) + if output.is_relative_to(root) or root.is_relative_to(output): + raise ValueError("Output and dataset must be disjoint directory trees") + if (output / "target_access.json").exists(): + raise ValueError( + "This run already started target scoring; use a separately identified pilot run" + ) + protocol = json.loads((output / "protocol.json").read_text()) + seal = json.loads((output / "prediction_seal.json").read_text()) + if ( + seal["protocol_id"] != object_digest(protocol) + or seal["implementation_sha256"] != implementation() + ): + raise ValueError("Protocol or implementation changed after prediction sealing") + if seal["source_fit_sha256"] != digest(output / "source_fit.json"): + raise ValueError("Source fit changed after sealing") + cases, inventory = audit_dataset(root, protocol) + if inventory["inventory_id"] != seal["inventory_id"]: + raise ValueError("Dataset changed after sealing") + for entry in seal["predictions"].values(): + path = (output / entry["artifact"]).resolve() + if ( + not path.is_relative_to((output / "private_predictions").resolve()) + or digest(path) != entry["sha256"] + ): + raise ValueError("Prediction artifact identity mismatch") + write_json( + output / "target_access.json", + { + "started_at": now(), + "prediction_seal_sha256": digest(output / "prediction_seal.json"), + "authorized_recordings": sorted(seal["predictions"]), + "purpose": "fixed public-data pilot scoring", + }, + ) + rows = [] + fallback_records = 0 + harmful_accepted_records = 0 + for case in (c for c in cases if c.motion == "twist"): + entry = seal["predictions"][case.path.name] + with np.load(output / entry["artifact"], allow_pickle=False) as arrays: + inputs = Inputs( + arrays["times"], + np.empty((0, case.markers, 3)), + np.empty((0, 2, 3)), + arrays["order"], + arrays["corners"], + int(arrays["cutoff"]), + float(arrays["times"][0]), + float(arrays["scale"]), + ) + truth = scoring_view(case, inputs) + case_scores = {} + for arm in ARMS: + mean, variance = arrays[f"{arm}_mean"], arrays[f"{arm}_variance"] + case_scores[arm] = score(mean, variance, truth, inputs) + rows.append( + { + "recording": case.path.name, + "specimen": case.specimen, + "material": case.material, + "speed": case.speed, + "grasp": case.grasp, + "arm": arm, + "guard_accepted": entry["guard_accepts"], + **case_scores[arm], + } + ) + if not entry["guard_accepts"]: + fallback_records += 1 + for field in ("mean", "variance"): + if not np.array_equal( + arrays[f"guarded_bayesian_physics_{field}"], + arrays[f"nominal_physics_{field}"], + ): + raise ValueError("Exact fallback violated") + if ( + case_scores["guarded_bayesian_physics"] + != case_scores["nominal_physics"] + ): + raise ValueError("Exact fallback score violated") + elif ( + case_scores["guarded_bayesian_physics"]["rmse_mm"] + > case_scores["nominal_physics"]["rmse_mm"] + ): + harmful_accepted_records += 1 + table, metrics = aggregate(rows, protocol) + metrics.update( + { + "fallback_recordings": fallback_records, + "accepted_recordings": 32 - fallback_records, + "harmful_accepted_recordings_vs_nominal": harmful_accepted_records, + "exact_fallback_violations": 0, + "target_recordings": 32, + "evidence_class": protocol["evidence_class"], + "paper_claim_authorized": False, + } + ) + save_csv(output / "target_scores.csv", rows) + save_csv(output / "specimen_scores.csv", table) + write_json(output / "metrics.json", metrics) + manifest = json.loads((output / "run_manifest.json").read_text()) + manifest.update( + { + "completed_at": now(), + "target_numeric_outcomes_read": True, + "prediction_seal_sha256": digest(output / "prediction_seal.json"), + "metrics_sha256": digest(output / "metrics.json"), + "status": "completed-pilot-not-claim-promoted", + } + ) + write_json(output / "run_manifest.json", manifest) + report = (output / "report.md").read_text() + report += "\n## Held-out twisting results\n\n" + report += "| Arm | Specimen-balanced RMSE [mm] | Coordinate NLL | 90% coverage | Full width [mm] |\n" + report += "| --- | ---: | ---: | ---: | ---: |\n" + for arm in ARMS: + values = metrics["arms"][arm] + report += ( + f"| {arm} | {values['rmse_mm']:.4f} | {values['coordinate_nll']:.4f} | " + f"{100 * values['coordinate_90_coverage']:.2f}% | {values['mean_full_90_width_mm']:.4f} |\n" + ) + report += ( + f"\nGuarded candidate: {32 - fallback_records}/32 accepted; {fallback_records}/32 exact fallbacks. " + f"{harmful_accepted_records} accepted records worsen RMSE versus nominal physics.\n\n" + "Intervals in metrics.json resample eight specimens, with a four-material sensitivity check. " + "They are exploratory, non-simultaneous intervals with few clusters. " + "Diagonal moment-matched Gaussian scores do not validate joint trajectory covariance.\n\n" + "The primary endpoint is free-marker Euclidean RMSE, not the paper's mass-matrix metric. " + "Known future measured corner positions are prescribed inputs; this is not command-conditioned " + "or fully online forecasting. No unseen-object, material-identification, causal or safety claim follows.\n" + ) + (output / "report.md").write_text(report) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--dataset-root", type=Path, required=True) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument( + "--stage", choices=("inventory", "source", "predict", "score"), default="source" + ) + parser.add_argument("--workers", type=int, default=4) + args = parser.parse_args() + if not 1 <= args.workers <= 8: + parser.error("workers must be between 1 and 8") + try: + if args.stage == "score": + score_run(args.dataset_root, args.output) + else: + protocol = json.loads((HERE / "protocol.json").read_text()) + prepare(args.dataset_root, args.output, protocol, args.stage, args.workers) + except Exception as exc: + if args.output.is_dir() and not args.output.resolve().is_relative_to( + args.dataset_root.resolve() + ): + write_json( + args.output / "failure.json", + { + "failed_at": now(), + "stage": args.stage, + "exception": type(exc).__name__, + "message": str(exc), + "target_scoring_started": ( + args.output / "target_access.json" + ).exists(), + "scientific_decision": "not-evaluated-or-incomplete; no claim", + }, + ) + traceback.print_exc() + return 1 + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/test_tracking_cloth_active_probe_run_v1.py b/tests/test_tracking_cloth_active_probe_run_v1.py new file mode 100644 index 000000000..ef37c529b --- /dev/null +++ b/tests/test_tracking_cloth_active_probe_run_v1.py @@ -0,0 +1,138 @@ +"""Synthetic integration contracts for the active-probe cloth runner.""" + +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np + +from experiments.tracking_cloth_deformation_v1.active_probe_run import ( + active_mask, + arm_specs, + array_digest, + belief_digest, + build_belief_arms, + calibrated_residuals, + loss_vector, + posterior_temperature, + validate_protocol, + weighted_belief, +) +from experiments.tracking_cloth_deformation_v1.data import Inputs +from experiments.tracking_cloth_deformation_v1.model import Predictions + +BASE = Path(__file__).resolve().parents[1] / "experiments/tracking_cloth_deformation_v1" + + +def protocol() -> dict: + value = json.loads((BASE / "active_probe_protocol.json").read_text()) + validate_protocol(value) + return value + + +def prediction() -> Predictions: + time_count, markers, models = 8, 4, 9 + times = np.arange(time_count, dtype=np.float64) + cutoff = 2 + prefix = np.zeros((cutoff + 1, markers, 3), dtype=np.float64) + boundary = np.zeros((time_count, 2, 3), dtype=np.float64) + inputs = Inputs( + times, + prefix, + boundary, + np.arange(markers), + np.array([0, markers - 1]), + cutoff, + 0.0, + 1.0, + ) + bank = np.zeros((models, time_count, markers, 3), dtype=np.float64) + for model in range(models): + bank[model, :, 1:-1, 0] = model * np.arange(time_count)[:, None] * 0.001 + return Predictions(inputs, bank[4].copy(), bank) + + +def test_protocol_registers_decisive_budget_and_information_boundary() -> None: + value = protocol() + assert value["fold_rule"] == "leave-one-material-out" + assert value["probe_budgets"] == [0, 1, 2, 4] + assert value["primary_budget"] == 1 + assert value["held_material_candidate_inputs_used_for_selection"] is False + assert value["held_material_twist_inputs_used_for_selection"] is False + assert value["paper_claim_authorized"] is False + + +def test_weighted_belief_and_content_digests() -> None: + candidate = prediction() + weights = np.ones(9) / 9 + mean, variance = weighted_belief(candidate, weights, [1e-4] * 3) + assert mean.shape == variance.shape == candidate.nominal.shape + assert np.all(variance > 0.0) + assert len(array_digest(mean)) == 64 + assert len(belief_digest(mean, variance)) == 64 + + +def test_complete_arm_roster_and_common_endpoint_parity() -> None: + value = protocol() + candidate = prediction() + weights = (np.ones(9) / 9).tolist() + states = { + policy: { + str(budget): { + "weights": weights, + "selected_actions": value["probe_conditions"][:budget], + "steps": [], + } + for budget in value["probe_budgets"] + } + for policy in value["probe_policies"] + } + specimen = { + "policy_states": states, + "single_probe_weights": { + condition: weights for condition in value["probe_conditions"] + }, + } + fold = { + "prior_weights": weights, + "source_residual_variance_m2": { + "bayesian": [1e-4] * 3, + "nominal_physics": [2e-4] * 3, + "last_residual": [3e-4] * 3, + }, + } + arms = build_belief_arms( + candidate, + prefix_last=candidate.inputs.prefix[-1], + boundary=candidate.inputs.boundary, + fold=fold, + specimen=specimen, + protocol=value, + ) + assert set(arms) == set(arm_specs(value)) + assert len(arms) == 18 + for budget in (0, 4): + reference = arms[f"fixed_order_k{budget}"] + for policy in ("parameter_information", "task_directed"): + np.testing.assert_array_equal(reference[0], arms[f"{policy}_k{budget}"][0]) + np.testing.assert_array_equal(reference[1], arms[f"{policy}_k{budget}"][1]) + + +def test_loss_temperature_residual_calibration_and_mask() -> None: + value = protocol() + candidate = prediction() + truth = candidate.bank[3].copy() + losses = loss_vector(candidate, truth) + assert losses[3] == 0.0 + assert posterior_temperature(np.stack([losses, losses]), 0.001) > 0.0 + residual = calibrated_residuals([(candidate, truth)], np.ones(9) / 9, value) + assert set(residual) == { + "bayesian", + "nominal_physics", + "last_residual", + } + assert all(number > 0.0 for values in residual.values() for number in values) + mask = active_mask(candidate.inputs) + assert mask.shape == (8, 4) + assert int(mask.sum()) == 10 diff --git a/tests/test_tracking_cloth_active_probe_v1.py b/tests/test_tracking_cloth_active_probe_v1.py new file mode 100644 index 000000000..9337c0740 --- /dev/null +++ b/tests/test_tracking_cloth_active_probe_v1.py @@ -0,0 +1,227 @@ +"""Synthetic contracts for outcome-blind active probe selection.""" + +from __future__ import annotations + +from collections.abc import Iterator, Mapping + +import numpy as np +import pytest + +from experiments.tracking_cloth_deformation_v1.active_probe import ( + entropy, + model_spread, + normalize_weights, + pairwise_trajectory_mse, + parameter_information_utility, + pseudo_posteriors, + select_action, + simulate_policy, + task_variance_reduction_utility, + update_weights, + weights_from_records, +) + + +class AccessLog(Mapping[str, np.ndarray]): + def __init__(self, values: dict[str, np.ndarray]): + self.values = values + self.accessed: list[str] = [] + + def __getitem__(self, key: str) -> np.ndarray: + self.accessed.append(key) + return self.values[key] + + def __iter__(self) -> Iterator[str]: + return iter(self.values) + + def __len__(self) -> int: + return len(self.values) + + +def distances(scale_a: float, scale_b: float = 0.0) -> np.ndarray: + # Three models; model 0 and 1 separate in the first mode, 1 and 2 in the second. + vectors = np.array([[0.0, 0.0], [scale_a, 0.0], [scale_a, scale_b]]) + delta = vectors[:, None] - vectors[None] + return np.sum(delta * delta, axis=-1) + + +def test_weights_and_updates_are_normalized_and_immutable() -> None: + weights = normalize_weights([1.0, 2.0, 0.0]) + assert weights.sum() == pytest.approx(1.0) + assert np.all(weights > 0.0) + with pytest.raises(ValueError): + weights[0] = 0.0 + updated = update_weights(weights, [0.0, 2.0, 4.0], 1.0) + assert updated[0] > weights[0] + fitted = weights_from_records([[0.0, 1.0, 2.0], [0.0, 1.0, 3.0]], 1.0) + assert fitted[0] > fitted[1] > fitted[2] + + +def test_pairwise_trajectory_metric_matches_point_squared_error() -> None: + bank = np.zeros((3, 2, 2, 3)) + bank[1, :, :, 0] = 2.0 + bank[2, :, :, 1] = 3.0 + valid = np.array([[False, True], [True, False]]) + matrix = pairwise_trajectory_mse(bank, valid) + assert matrix[0, 1] == pytest.approx(4.0) + assert matrix[0, 2] == pytest.approx(9.0) + assert matrix[1, 2] == pytest.approx(13.0) + np.testing.assert_array_equal(np.diag(matrix), 0.0) + + +def test_information_and_task_objectives_can_choose_different_probes() -> None: + weights = np.array([0.45, 0.45, 0.10]) + # Probe a separates the high-mass pair; probe b isolates the low-mass model. + probe_a = distances(3.0, 0.0) + probe_b = distances(0.0, 8.0) + # The held-out task only distinguishes model 2 from the first two. + target = distances(0.0, 20.0) + parameter = { + action: parameter_information_utility(weights, matrix, 1.0) + for action, matrix in {"a": probe_a, "b": probe_b}.items() + } + task = { + action: task_variance_reduction_utility(weights, matrix, target, 1.0) + for action, matrix in {"a": probe_a, "b": probe_b}.items() + } + # Verify the objectives are genuinely different rather than hard-code an order. + assert parameter["a"] > parameter["b"] + assert task["b"] > task["a"] + assert ( + select_action( + policy="parameter_information", + weights=weights, + remaining_actions=("a", "b"), + probe_distances={"a": probe_a, "b": probe_b}, + target_distance=target, + temperature=1.0, + fixed_order=("a", "b"), + )[0] + == "a" + ) + assert ( + select_action( + policy="task_directed", + weights=weights, + remaining_actions=("a", "b"), + probe_distances={"a": probe_a, "b": probe_b}, + target_distance=target, + temperature=1.0, + fixed_order=("a", "b"), + )[0] + == "b" + ) + + +def test_task_utility_is_expected_target_spread_contraction() -> None: + weights = normalize_weights([0.4, 0.4, 0.2]) + probe = distances(0.0, 5.0) + target = distances(0.0, 10.0) + posteriors = pseudo_posteriors(weights, probe, 0.5) + expected = sum( + weights[index] * model_spread(posteriors[index], target) for index in range(3) + ) + manual = 1.0 - expected / model_spread(weights, target) + assert task_variance_reduction_utility( + weights, probe, target, 0.5 + ) == pytest.approx(manual) + assert entropy(weights) > 0.0 + + +def test_policy_consumes_only_selected_outcomes() -> None: + probe = { + "a": distances(3.0), + "b": distances(0.0, 8.0), + "c": distances(1.0, 1.0), + } + losses = AccessLog( + { + "a": np.array([0.0, 1.0, 2.0]), + "b": np.array([2.0, 1.0, 0.0]), + "c": np.array([1.0, 0.0, 1.0]), + } + ) + states = simulate_policy( + policy="task_directed", + initial_weights=[0.45, 0.45, 0.10], + probe_distances=probe, + target_distance=distances(0.0, 20.0), + observed_losses=losses, + temperature=1.0, + fixed_order=("a", "b", "c"), + budgets=(0, 1, 2), + ) + assert losses.accessed == list(states[2].selected_actions) + assert len(losses.accessed) == 2 + assert states[0].selected_actions == () + assert len(states[1].selected_actions) == 1 + assert len(states[2].selected_actions) == 2 + + +def test_all_probe_posterior_is_order_invariant() -> None: + probe = { + "a": distances(3.0), + "b": distances(0.0, 8.0), + "c": distances(1.0, 1.0), + } + values = { + "a": np.array([0.0, 1.0, 2.0]), + "b": np.array([2.0, 1.0, 0.0]), + "c": np.array([1.0, 0.0, 1.0]), + } + finals = [] + for policy in ("fixed_order", "parameter_information", "task_directed"): + states = simulate_policy( + policy=policy, + initial_weights=[0.4, 0.4, 0.2], + probe_distances=probe, + target_distance=distances(0.0, 20.0), + observed_losses=values, + temperature=1.0, + fixed_order=("c", "a", "b"), + budgets=(0, 1, 2, 3), + ) + finals.append(states[3].weights) + np.testing.assert_allclose(finals[0], finals[1], atol=1e-15) + np.testing.assert_allclose(finals[0], finals[2], atol=1e-15) + + +def test_deterministic_tie_break_and_fixed_order() -> None: + zeros = np.zeros((3, 3)) + distances_by_action = {"z": zeros, "a": zeros} + selected, _ = select_action( + policy="task_directed", + weights=[1, 1, 1], + remaining_actions=("z", "a"), + probe_distances=distances_by_action, + target_distance=zeros, + temperature=1.0, + fixed_order=("z", "a"), + ) + assert selected == "a" + selected, utilities = select_action( + policy="fixed_order", + weights=[1, 1, 1], + remaining_actions=("a", "z"), + probe_distances=distances_by_action, + target_distance=zeros, + temperature=1.0, + fixed_order=("z", "a"), + ) + assert selected == "z" + assert utilities == {"a": None, "z": None} + + +@pytest.mark.parametrize( + "call", + [ + lambda: normalize_weights([0.0, 0.0]), + lambda: normalize_weights([1.0]), + lambda: update_weights([1, 1], [0, -1], 1.0), + lambda: pseudo_posteriors([1, 1], [[0, 1], [2, 0]], 1.0), + lambda: pairwise_trajectory_mse(np.zeros((2, 1, 1, 2)), [[True]]), + ], +) +def test_invalid_inputs_fail_closed(call) -> None: + with pytest.raises(ValueError): + call() diff --git a/tests/test_tracking_cloth_deformation_v1.py b/tests/test_tracking_cloth_deformation_v1.py new file mode 100644 index 000000000..ce393405b --- /dev/null +++ b/tests/test_tracking_cloth_deformation_v1.py @@ -0,0 +1,312 @@ +"""Synthetic-only data, leakage, numerical and evaluation contracts.""" + +from __future__ import annotations + +import csv +import hashlib +import io +import itertools +import json +import zipfile +from pathlib import Path + +import numpy as np +import pytest + +from experiments.tracking_cloth_deformation_v1.data import ( + Case, + audit_dataset, + input_view, + layout, + scoring_view, +) +from experiments.tracking_cloth_deformation_v1.model import ( + ARMS, + Predictions, + complete_beliefs, + masks, + parameter_bank, + predict, + score, + source_weights, +) +from experiments.tracking_cloth_deformation_v1.run import aggregate, prepare, score_run + +BASE = Path(__file__).resolve().parents[1] / "experiments/tracking_cloth_deformation_v1" + + +def config(): + protocol = json.loads((BASE / "protocol.json").read_text()) + protocol.update( + { + "prefix_seconds": 0.2, + "forecast_seconds": 0.2, + "stiffness_per_mass": [400.0], + "damping_per_mass": [2.0], + "integration_substeps": 2, + "bootstrap_repetitions": 100, + } + ) + return protocol + + +def initial_grid(size="A3"): + rows, cols = (5, 4) if size == "A2" else (4, 3) + spacing = 0.12 + return np.array( + [ + [spacing * c, 0.0, 1.0 - spacing * r] + for r in range(rows) + for c in range(cols) + ] + ) + + +def csv_text(size="A3", poison_future=False, missing=False): + first = initial_grid(size) + stream = io.StringIO() + writer = csv.writer(stream) + writer.writerow(["Format Version", "1.0", "Length Units", "Meters"]) + writer.writerow(["", "", "Position"]) + writer.writerow(["", "", "Marker IDs"]) + writer.writerow(["", "", "Position"]) + writer.writerow(["Frame", "Time", *(["X", "Y", "Z"] * len(first))]) + for i in range(61): + t = i / 120 + positions = first.copy() + positions[:, 1] += 0.002 * np.sin(2 * np.pi * t) + values = positions.reshape(-1).astype(object) + if missing and i == 40: + values[3 * 5 : 3 * 5 + 3] = "" + if poison_future and t > 0.2 + 1e-8: + # Only non-driven markers are poisoned, proving no numeric conversion. + for marker in range(len(first)): + if marker not in (0, 2 if size == "A3" else 3): + values[3 * marker : 3 * marker + 3] = "UNOPENED_TARGET" + writer.writerow([i, f"{t:.9f}", *values]) + return stream.getvalue() + + +def case_file(tmp_path, text=None, motion="twist"): + path = tmp_path / f"cotton_A3_{motion}_fast_hands.csv" + path.write_text(csv_text() if text is None else text) + return Case(path, "cotton", "A3", motion, "fast", "hands") + + +@pytest.fixture +def dataset(tmp_path): + root = tmp_path / "dataset" + root.mkdir() + payloads = {} + for material, size, motion, speed, grasp in itertools.product( + ("cotton", "denim", "polyester", "wool"), + ("A2", "A3"), + ("shake", "twist"), + ("fast", "slow"), + ("hands", "hanger"), + ): + name = f"Free-hanging/{material}_{size}_{motion}_{speed}_{grasp}.csv" + payloads[name] = csv_text(size) + for i in range(56): + payloads[f"Reserved/unused_{i}.csv"] = "reserved; never numerically parsed\n" + for folder in ("Free-hanging", "Tablecloth", "Hitting", "Self-collision"): + payloads[f"{folder}/read_data.m"] = "% Synthetic fixture, no dataset content.\n" + payloads["License.txt"] = "Synthetic fixture: CC BY-NC-SA 4.0\n" + for name, content in payloads.items(): + path = root / name + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(content) + archive = root / "dataset.zip" + with zipfile.ZipFile(archive, "w", zipfile.ZIP_DEFLATED) as zipped: + for name, content in payloads.items(): + zipped.writestr("dataset/" + name, content) + protocol = config() + protocol["archive_md5"] = hashlib.md5( + archive.read_bytes(), usedforsecurity=False + ).hexdigest() + return root, protocol + + +def test_complete_roster_inventory_without_numeric_target_reads(dataset): + root, protocol = dataset + cases, inventory = audit_dataset(root, protocol) + assert len(cases) == 64 + assert inventory["source_count"] == inventory["target_count"] == 32 + assert inventory["unused_count"] == 56 + assert inventory["target_numeric_outcomes_read"] is False + + +def test_tampered_extraction_rejected(dataset): + root, protocol = dataset + next(root.rglob("*twist*.csv")).write_text("tampered") + with pytest.raises(ValueError, match="Extracted bytes"): + audit_dataset(root, protocol) + + +def test_duplicate_basename_rejected(dataset): + root, protocol = dataset + duplicate = root / "dup" + duplicate.mkdir() + original = next(root.rglob("*twist*.csv")) + (duplicate / original.name).write_bytes(original.read_bytes()) + protocol["csv_count"] += 1 + with pytest.raises(ValueError, match="duplicate"): + audit_dataset(root, protocol) + + +def test_wrong_archive_rejected(dataset): + root, protocol = dataset + protocol["archive_md5"] = "0" * 32 + with pytest.raises(ValueError, match="published MD5"): + audit_dataset(root, protocol) + + +def test_canonical_grid_and_driven_corners(): + permutation = np.random.default_rng(8).permutation(20) + order, corners = layout(initial_grid("A2")[permutation], "A2") + np.testing.assert_array_equal(permutation[order], np.arange(20)) + np.testing.assert_array_equal(corners, [0, 3]) + + +def test_future_free_marker_values_never_enter_prediction(tmp_path): + case = case_file(tmp_path) + protocol = config() + clean = input_view(case, protocol, 1.0) + before = predict(clean, protocol) + case.path.write_text(csv_text(poison_future=True)) + closed = input_view(case, protocol, 1.0) + after = predict(closed, protocol) + np.testing.assert_array_equal(before.nominal, after.nominal) + np.testing.assert_array_equal(before.bank, after.bank) + with pytest.raises(ValueError, match="UNOPENED_TARGET"): + scoring_view(case, closed) + + +def test_causal_missing_values_not_filled_from_future(tmp_path): + case = case_file(tmp_path, csv_text(missing=True)) + inputs = input_view(case, config(), 1.0) + truth = scoring_view(case, inputs) + assert np.isnan(truth).any() + assert np.isfinite(inputs.prefix).all() + valid = masks(inputs, truth) + assert not np.any(valid[:, inputs.corners]) + assert not np.any(valid[: inputs.cutoff + 1]) + + +def test_bad_timestamps_rejected(tmp_path): + case = case_file(tmp_path, csv_text().replace("1,0.008333333", "1,0.000000000")) + with pytest.raises(ValueError, match="order"): + input_view(case, config(), 1.0) + + +def test_gibbs_weights_prefer_lower_source_loss(): + weights = source_weights(np.array([[0.001, 0.1], [0.002, 0.2]]), 0.001) + assert weights[0] > 0.999 + assert weights.sum() == pytest.approx(1.0) + + +def test_exact_fallback_includes_covariance(tmp_path): + inputs = input_view(case_file(tmp_path), config(), 1.0) + prediction = predict(inputs, config()) + fit = { + "source_posterior_weights": [1.0], + "guard_accepts": False, + "source_residual_variance_m2": {arm: [1e-5, 2e-5, 3e-5] for arm in ARMS[:-1]}, + } + beliefs = complete_beliefs(prediction, fit, config()) + assert beliefs["guarded_bayesian_physics"] is beliefs["nominal_physics"] + np.testing.assert_array_equal( + prediction.nominal[:, inputs.corners], inputs.boundary + ) + np.testing.assert_array_equal( + prediction.bank[0][:, inputs.corners], inputs.boundary + ) + + +def test_coordinate_score_ignores_driven_markers(tmp_path): + inputs = input_view(case_file(tmp_path), config(), 1.0) + truth = scoring_view(case_file(tmp_path), inputs) + mean = truth.copy() + mean[:, inputs.corners] += 10000 + values = score(mean, np.full_like(mean, 1e-4), truth, inputs) + assert values["rmse_mm"] == 0 + assert values["coordinate_90_coverage"] == 1 + + +def test_posterior_total_variance_includes_parameter_spread(tmp_path): + protocol = config() + protocol["stiffness_per_mass"] = [100.0, 400.0] + inputs = input_view(case_file(tmp_path), protocol, 1.0) + nominal = np.zeros((len(inputs.times), 12, 3)) + prediction = Predictions(inputs, nominal, np.stack([nominal, nominal + 2])) + fit = { + "source_posterior_weights": [0.5, 0.5], + "guard_accepts": True, + "source_residual_variance_m2": {arm: [1.0] * 3 for arm in ARMS[:-1]}, + } + beliefs = complete_beliefs(prediction, fit, protocol) + assert np.all(beliefs["bayesian_physics"][0] == 1) + assert np.all(beliefs["bayesian_physics"][1] == 2) + assert beliefs["guarded_bayesian_physics"] is beliefs["bayesian_physics"] + assert parameter_bank(protocol) == [(100.0, 2.0), (400.0, 2.0)] + + +def test_source_only_run_does_not_touch_target_numeric_payload(dataset, tmp_path): + root, protocol = dataset + output = tmp_path / "source-output" + prepare(root, output, protocol, "source") + assert (output / "source_fit.json").is_file() + assert not (output / "private_predictions").exists() + assert not (output / "target_access.json").exists() + assert not (output / "metrics.json").exists() + assert ( + json.loads((output / "source_fit.json").read_text())["target_outcomes_used"] + is False + ) + + +def test_complete_predict_seal_score_cycle(dataset, tmp_path): + root, protocol = dataset + output = tmp_path / "output" + prepare(root, output, protocol, "predict") + assert (output / "prediction_seal.json").is_file() + assert not (output / "target_access.json").exists() + assert not (output / "metrics.json").exists() + score_run(root, output) + metrics = json.loads((output / "metrics.json").read_text()) + assert metrics["target_recordings"] == 32 + assert metrics["paper_claim_authorized"] is False + assert metrics["exact_fallback_violations"] == 0 + assert set(metrics["arms"]) == set(ARMS) + with pytest.raises(ValueError, match="already started"): + score_run(root, output) + + +def test_output_cannot_be_dataset_descendant(dataset): + root, protocol = dataset + with pytest.raises(ValueError, match="disjoint"): + prepare(root, root / "outputs", protocol, "inventory") + + +def test_modified_source_fit_stops_target_opening(dataset, tmp_path): + root, protocol = dataset + output = tmp_path / "altered-output" + prepare(root, output, protocol, "predict") + (output / "source_fit.json").write_text("{}") + with pytest.raises(ValueError, match="Source fit changed"): + score_run(root, output) + assert not (output / "target_access.json").exists() + + +def test_incomplete_result_never_pooled(): + with pytest.raises(ValueError, match="Incomplete roster"): + aggregate([], config()) + + +def test_checked_in_protocol_matches_requested_release(): + protocol = json.loads((BASE / "protocol.json").read_text()) + assert protocol["dataset_record"] == "14644526" + assert protocol["archive_md5"] == "b4868b702f8a42b2ea1069d0f1a3b8f6" + assert protocol["raw_data_upload"] is False + assert protocol["paper_claim_authorized"] is False + assert "NC" in protocol["license_policy"]