|
| 1 | +#!/usr/bin/env python3 |
| 2 | +"""Physical CuPy/Torch validation for Panel P1 Stage A (Issue #93 / PR #119). |
| 3 | +
|
| 4 | +This is a correctness/backend acceptance script, not a performance benchmark. |
| 5 | +It runs the behavior-preserving Stage-A panel refactor on deterministic panel |
| 6 | +data, compares CuPy/Torch results with the NumPy reference, and records exact |
| 7 | +source/environment provenance. |
| 8 | +""" |
| 9 | + |
| 10 | +from __future__ import annotations |
| 11 | + |
| 12 | +import argparse |
| 13 | +import importlib.metadata |
| 14 | +import json |
| 15 | +import platform |
| 16 | +import subprocess |
| 17 | +from datetime import datetime, timezone |
| 18 | +from pathlib import Path |
| 19 | + |
| 20 | +import numpy as np |
| 21 | + |
| 22 | +from statgpu.backends import _to_numpy |
| 23 | +from statgpu.panel import ( |
| 24 | + BetweenOLS, |
| 25 | + FamaMacBeth, |
| 26 | + FirstDifferenceOLS, |
| 27 | + PanelOLS, |
| 28 | + PooledOLS, |
| 29 | + RandomEffects, |
| 30 | +) |
| 31 | + |
| 32 | + |
| 33 | +def _git_sha() -> str: |
| 34 | + try: |
| 35 | + return subprocess.check_output( |
| 36 | + ["git", "rev-parse", "HEAD"], text=True |
| 37 | + ).strip() |
| 38 | + except (OSError, subprocess.CalledProcessError): |
| 39 | + return "unknown" |
| 40 | + |
| 41 | + |
| 42 | +def _git_status_porcelain() -> str: |
| 43 | + try: |
| 44 | + return subprocess.check_output( |
| 45 | + ["git", "status", "--porcelain"], text=True |
| 46 | + ) |
| 47 | + except (OSError, subprocess.CalledProcessError) as exc: |
| 48 | + raise RuntimeError("unable to verify working-tree cleanliness") from exc |
| 49 | + |
| 50 | + |
| 51 | +def _version(name: str): |
| 52 | + try: |
| 53 | + return importlib.metadata.version(name) |
| 54 | + except importlib.metadata.PackageNotFoundError: |
| 55 | + return None |
| 56 | + |
| 57 | + |
| 58 | +def _dataset(): |
| 59 | + rng = np.random.default_rng(20260807) |
| 60 | + n_entities, n_times = 8, 6 |
| 61 | + entity = np.repeat(np.arange(n_entities), n_times) |
| 62 | + time = np.tile(np.arange(n_times), n_entities) |
| 63 | + X = rng.normal(size=(entity.size, 2)) |
| 64 | + entity_effect = np.repeat(rng.normal(scale=0.6, size=n_entities), n_times) |
| 65 | + time_effect = np.tile(np.linspace(-0.25, 0.30, n_times), n_entities) |
| 66 | + y = ( |
| 67 | + 0.4 |
| 68 | + + 1.15 * X[:, 0] |
| 69 | + - 0.65 * X[:, 1] |
| 70 | + + entity_effect |
| 71 | + + 0.20 * time_effect |
| 72 | + + rng.normal(scale=0.16, size=entity.size) |
| 73 | + ) |
| 74 | + return X.astype(np.float64), y.astype(np.float64), entity, time |
| 75 | + |
| 76 | + |
| 77 | +def _to_backend(X, y, backend): |
| 78 | + if backend == "numpy": |
| 79 | + return X, y |
| 80 | + if backend == "cupy": |
| 81 | + import cupy as cp |
| 82 | + |
| 83 | + return cp.asarray(X), cp.asarray(y) |
| 84 | + if backend == "torch": |
| 85 | + import torch |
| 86 | + |
| 87 | + return ( |
| 88 | + torch.as_tensor(X, dtype=torch.float64, device="cuda"), |
| 89 | + torch.as_tensor(y, dtype=torch.float64, device="cuda"), |
| 90 | + ) |
| 91 | + raise ValueError(backend) |
| 92 | + |
| 93 | + |
| 94 | +def _device_arg(backend): |
| 95 | + return {"numpy": "cpu", "cupy": "cuda", "torch": "torch"}[backend] |
| 96 | + |
| 97 | + |
| 98 | +def _backend_name(model): |
| 99 | + if isinstance(model, FamaMacBeth): |
| 100 | + return model._backend_name |
| 101 | + return model._get_backend(backend="auto").name |
| 102 | + |
| 103 | + |
| 104 | +def _array(value): |
| 105 | + return np.asarray(_to_numpy(value), dtype=np.float64) |
| 106 | + |
| 107 | + |
| 108 | +def _snapshot(model, prediction): |
| 109 | + payload = { |
| 110 | + "coef": _array(model.coef_).ravel(), |
| 111 | + "bse": _array(model.bse_).ravel(), |
| 112 | + "tvalues": _array(model.tvalues_).ravel(), |
| 113 | + "pvalues": _array(model.pvalues_).ravel(), |
| 114 | + "conf_int": _array(model.conf_int_), |
| 115 | + "prediction": _array(prediction).ravel(), |
| 116 | + "df_resid": int(model.df_resid), |
| 117 | + "nobs": int(model.nobs), |
| 118 | + } |
| 119 | + if hasattr(model, "rsquared"): |
| 120 | + payload["rsquared"] = float(model.rsquared) |
| 121 | + if hasattr(model, "rsquared_within") and model.rsquared_within is not None: |
| 122 | + payload["rsquared_within"] = float(model.rsquared_within) |
| 123 | + if hasattr(model, "theta_") and model.theta_ is not None: |
| 124 | + payload["theta"] = float(model.theta_) |
| 125 | + if hasattr(model, "variance_components_") and model.variance_components_ is not None: |
| 126 | + payload["variance_components"] = { |
| 127 | + key: float(value) |
| 128 | + for key, value in model.variance_components_.items() |
| 129 | + } |
| 130 | + if hasattr(model, "n_periods"): |
| 131 | + payload["n_periods"] = int(model.n_periods) |
| 132 | + return payload |
| 133 | + |
| 134 | + |
| 135 | +def _cases(X, y, entity, time, backend): |
| 136 | + Xb, yb = _to_backend(X, y, backend) |
| 137 | + device = _device_arg(backend) |
| 138 | + two_way = np.column_stack([entity, time]) |
| 139 | + |
| 140 | + cases = {} |
| 141 | + |
| 142 | + model = PooledOLS(cov_type="nonrobust", device=device).fit(Xb, yb) |
| 143 | + cases["pooled_nonrobust"] = (model, model.predict(X[:5])) |
| 144 | + |
| 145 | + model = PooledOLS(cov_type="robust", device=device).fit(Xb, yb) |
| 146 | + cases["pooled_robust"] = (model, model.predict(X[:5])) |
| 147 | + |
| 148 | + model = PooledOLS(cov_type="clustered", device=device).fit( |
| 149 | + Xb, yb, cluster=entity |
| 150 | + ) |
| 151 | + cases["pooled_clustered"] = (model, model.predict(X[:5])) |
| 152 | + |
| 153 | + model = PooledOLS(cov_type="hac", bandwidth=2, device=device).fit( |
| 154 | + Xb, yb, time_index=time |
| 155 | + ) |
| 156 | + cases["pooled_hac"] = (model, model.predict(X[:5])) |
| 157 | + |
| 158 | + model = BetweenOLS(cov_type="robust", device=device).fit( |
| 159 | + Xb, yb, entity_ids=entity |
| 160 | + ) |
| 161 | + cases["between_robust"] = (model, model.predict(X[:5])) |
| 162 | + |
| 163 | + model = FirstDifferenceOLS(cov_type="robust", device=device).fit( |
| 164 | + Xb, yb, entity_ids=entity, time_ids=time |
| 165 | + ) |
| 166 | + cases["first_difference_robust"] = (model, model.predict(X[:5])) |
| 167 | + |
| 168 | + model = PanelOLS(entity_effects=True, cov_type="robust", device=device).fit( |
| 169 | + Xb, yb, entity_ids=entity |
| 170 | + ) |
| 171 | + cases["panel_entity_robust"] = ( |
| 172 | + model, |
| 173 | + model.predict(X[:5], entity_ids=entity[:5]), |
| 174 | + ) |
| 175 | + |
| 176 | + model = PanelOLS( |
| 177 | + entity_effects=True, |
| 178 | + time_effects=True, |
| 179 | + cov_type="clustered", |
| 180 | + device=device, |
| 181 | + ).fit( |
| 182 | + Xb, |
| 183 | + yb, |
| 184 | + entity_ids=entity, |
| 185 | + time_ids=time, |
| 186 | + cluster=two_way, |
| 187 | + ) |
| 188 | + cases["panel_two_way_clustered"] = ( |
| 189 | + model, |
| 190 | + model.predict(X[:5], entity_ids=entity[:5], time_ids=time[:5]), |
| 191 | + ) |
| 192 | + |
| 193 | + model = RandomEffects(device=device).fit(Xb, yb, entity_ids=entity) |
| 194 | + cases["random_effects"] = (model, model.predict(X[:5])) |
| 195 | + |
| 196 | + model = FamaMacBeth(cov_type="newey-west", bandwidth=2, device=device).fit( |
| 197 | + Xb, yb, time_ids=time |
| 198 | + ) |
| 199 | + cases["fama_macbeth_newey_west"] = (model, model.predict(X[:5])) |
| 200 | + |
| 201 | + return cases |
| 202 | + |
| 203 | + |
| 204 | +def _compare(reference, candidate, *, rtol, atol): |
| 205 | + differences = {} |
| 206 | + for key in ("coef", "bse", "tvalues", "pvalues", "conf_int", "prediction"): |
| 207 | + np.testing.assert_allclose( |
| 208 | + candidate[key], reference[key], rtol=rtol, atol=atol, |
| 209 | + err_msg=f"mismatch in {key}", |
| 210 | + ) |
| 211 | + differences[key] = float( |
| 212 | + np.max(np.abs(np.asarray(candidate[key]) - np.asarray(reference[key]))) |
| 213 | + ) |
| 214 | + for key in ("df_resid", "nobs", "n_periods"): |
| 215 | + if key in reference: |
| 216 | + assert candidate[key] == reference[key], (key, candidate[key], reference[key]) |
| 217 | + for key in ("rsquared", "rsquared_within", "theta"): |
| 218 | + if key in reference: |
| 219 | + np.testing.assert_allclose(candidate[key], reference[key], rtol=rtol, atol=atol) |
| 220 | + differences[key] = float(abs(candidate[key] - reference[key])) |
| 221 | + if "variance_components" in reference: |
| 222 | + for name, value in reference["variance_components"].items(): |
| 223 | + np.testing.assert_allclose( |
| 224 | + candidate["variance_components"][name], value, rtol=rtol, atol=atol |
| 225 | + ) |
| 226 | + differences[f"variance_components.{name}"] = float( |
| 227 | + abs(candidate["variance_components"][name] - value) |
| 228 | + ) |
| 229 | + return differences |
| 230 | + |
| 231 | + |
| 232 | +def _environment(backends): |
| 233 | + gpu = None |
| 234 | + if "torch" in backends: |
| 235 | + import torch |
| 236 | + if not torch.cuda.is_available(): |
| 237 | + raise RuntimeError("Torch backend requested but CUDA is unavailable") |
| 238 | + gpu = torch.cuda.get_device_name(0) |
| 239 | + elif "cupy" in backends: |
| 240 | + import cupy as cp |
| 241 | + if cp.cuda.runtime.getDeviceCount() < 1: |
| 242 | + raise RuntimeError("CuPy backend requested but CUDA is unavailable") |
| 243 | + props = cp.cuda.runtime.getDeviceProperties(0) |
| 244 | + gpu = props["name"].decode() if isinstance(props["name"], bytes) else props["name"] |
| 245 | + return { |
| 246 | + "python": platform.python_version(), |
| 247 | + "platform": platform.platform(), |
| 248 | + "gpu": gpu, |
| 249 | + "packages": { |
| 250 | + name: _version(name) |
| 251 | + for name in ("statgpu", "numpy", "scipy", "cupy", "torch") |
| 252 | + }, |
| 253 | + } |
| 254 | + |
| 255 | + |
| 256 | +def main(): |
| 257 | + parser = argparse.ArgumentParser(description=__doc__) |
| 258 | + parser.add_argument("--out", type=Path, required=True) |
| 259 | + parser.add_argument( |
| 260 | + "--backends", default="cupy,torch", |
| 261 | + help="Comma-separated physical GPU backends to validate against NumPy.", |
| 262 | + ) |
| 263 | + parser.add_argument("--expected-sha", default=None) |
| 264 | + parser.add_argument("--rtol", type=float, default=5e-6) |
| 265 | + parser.add_argument("--atol", type=float, default=5e-7) |
| 266 | + args = parser.parse_args() |
| 267 | + |
| 268 | + backends = [value.strip() for value in args.backends.split(",") if value.strip()] |
| 269 | + if not backends or any(value not in {"cupy", "torch"} for value in backends): |
| 270 | + raise ValueError("--backends must contain cupy and/or torch") |
| 271 | + |
| 272 | + sha = _git_sha() |
| 273 | + if args.expected_sha is not None and sha != args.expected_sha: |
| 274 | + raise RuntimeError(f"wrong source head: {sha} != {args.expected_sha}") |
| 275 | + dirty = _git_status_porcelain() |
| 276 | + if dirty.strip(): |
| 277 | + raise RuntimeError( |
| 278 | + "physical acceptance requires a clean working tree; uncommitted changes:\n" |
| 279 | + + dirty |
| 280 | + ) |
| 281 | + |
| 282 | + X, y, entity, time = _dataset() |
| 283 | + reference_models = _cases(X, y, entity, time, "numpy") |
| 284 | + references = { |
| 285 | + name: _snapshot(model, prediction) |
| 286 | + for name, (model, prediction) in reference_models.items() |
| 287 | + } |
| 288 | + |
| 289 | + results = {} |
| 290 | + for backend in backends: |
| 291 | + backend_models = _cases(X, y, entity, time, backend) |
| 292 | + backend_results = {} |
| 293 | + for name, (model, prediction) in backend_models.items(): |
| 294 | + actual_backend = _backend_name(model) |
| 295 | + if actual_backend != backend: |
| 296 | + raise AssertionError( |
| 297 | + f"{name}: requested {backend}, executed {actual_backend}" |
| 298 | + ) |
| 299 | + snapshot = _snapshot(model, prediction) |
| 300 | + differences = _compare( |
| 301 | + references[name], snapshot, rtol=args.rtol, atol=args.atol |
| 302 | + ) |
| 303 | + backend_results[name] = { |
| 304 | + "status": "success", |
| 305 | + "executed_backend": actual_backend, |
| 306 | + "max_abs_differences": differences, |
| 307 | + } |
| 308 | + results[backend] = backend_results |
| 309 | + |
| 310 | + payload = { |
| 311 | + "schema_version": 1, |
| 312 | + "generated_at": datetime.now(timezone.utc).isoformat().replace("+00:00", "Z"), |
| 313 | + "git_sha": sha, |
| 314 | + "working_tree_clean": True, |
| 315 | + "status": "success", |
| 316 | + "environment": _environment(backends), |
| 317 | + "tolerances": {"rtol": args.rtol, "atol": args.atol}, |
| 318 | + "backends": results, |
| 319 | + } |
| 320 | + args.out.parent.mkdir(parents=True, exist_ok=True) |
| 321 | + args.out.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8") |
| 322 | + print(json.dumps(payload, indent=2)) |
| 323 | + print(f"PASS — Panel Stage A physical GPU validation: {args.out}") |
| 324 | + |
| 325 | + |
| 326 | +if __name__ == "__main__": |
| 327 | + main() |
0 commit comments