Add AGENTS paper registry and backend-selectable estimation path (CPU/CUDA/MLX/OpenCL)

README.md

@@ -69,6 +69,8 @@ print(process_dimensions.best)
 - Standalone HTML reports for saved fit or dimensionality diagnostics.
 - CLI commands for simulation and fitting.
 - Rust core crate with the same likelihood and gradient formulas.
+- Optional fit-time compute backend selection (`cpu`, `cuda`, `mlx`, `opencl`)
+  for model estimation.
 
 ## Install
 
@@ -101,6 +103,7 @@ fast-mlsirm fit \
   --responses runs/sim_001/responses.npy \
   --factors runs/sim_001/item_factor.csv \
   --model MLS2PLM \
+  --device cpu \
   --latent-dim 2 \
   --optimizer adam_lbfgs \
   --max-iter 100 \

docs/agents_papers.md

@@ -0,0 +1,16 @@
+# AGENTS.md referenced papers
+
+This repository task stores the paper set listed in `/home/runner/work/fast-mlsirm/fast-mlsirm/AGENTS.md`
+for traceable model/diagnostics work.
+
+1. Kang, I., & Jeon, M. (2025). *Multidimensional Latent Space Item Response Models: A Note on the Relativity of Conditional Dependence.* Psychometrika, 90(2), 799-826. https://doi.org/10.1017/psy.2025.5
+2. Jeon, M., Jin, I. H., Schweinberger, M., & Baugh, S. (2021). *Mapping Unobserved Item-Respondent Interactions: A Latent Space Item Response Model with Interaction Map.* Psychometrika, 86(2), 378-403. https://doi.org/10.1007/s11336-021-09762-5
+3. Molenaar, D., & Jeon, M. (2026). *Regularized Joint Maximum Likelihood Estimation of Latent Space Item Response Models.* Psychometrika, 91, 335-359. https://doi.org/10.1017/psy.2025.10068
+4. Tay, L., Ali, U. S., Drasgow, F., & Williams, B. (2011). *Fitting IRT Models to Dichotomous and Polytomous Data: Assessing the Relative Model-Data Fit of Ideal Point and Dominance Models.* Applied Psychological Measurement, 35(4), 280-295. https://doi.org/10.1177/0146621610390674
+5. Roberts, J. S., Donoghue, J. R., & Laughlin, J. E. (1998). *The Generalized Graded Unfolding Model: A General Parametric Item Response Model for Unfolding Graded Responses.* ETS Research Report Series. https://doi.org/10.1002/j.2333-8504.1998.tb01781.x
+6. Orlando, M., & Thissen, D. (2000). *Likelihood-Based Item-Fit Indices for Dichotomous Item Response Theory Models.* Applied Psychological Measurement, 24, 50-64.
+7. Maydeu-Olivares, A., & Joe, H. (2005). *Limited- and Full-Information Estimation and Goodness-of-Fit Testing in 2^n Contingency Tables.* Journal of the American Statistical Association, 100(471), 1009-1020. https://doi.org/10.1198/016214504000002069
+8. Drasgow, F., Levine, M. V., & Williams, E. A. (1985). *Appropriateness Measurement with Polychotomous Item Response Models and Standardized Indices.* British Journal of Mathematical and Statistical Psychology, 38(1), 67-86. https://doi.org/10.1111/j.2044-8317.1985.tb00817.x
+9. Fox, J.-P., & Glas, C. A. W. (2001). *Bayesian Estimation of a Multilevel IRT Model.* Psychometrika, 66, 271-288. https://doi.org/10.1007/BF02294839
+10. Bock, R. D., & Zimowski, M. F. (1997). *Multiple Group IRT.* In W. J. van der Linden & R. K. Hambleton (Eds.), *Handbook of Modern Item Response Theory.*
+11. Chalmers, R. P. (2012). *mirt: A Multidimensional Item Response Theory Package for the R Environment.* Journal of Statistical Software, 48(6). https://doi.org/10.18637/jss.v048.i06

pyproject.toml

@@ -14,6 +14,7 @@ dependencies = ["numpy>=1.24"]
 
 [project.optional-dependencies]
 dev = ["pytest>=8"]
+gpu = ["cupy-cuda12x>=13; platform_system == 'Linux'", "mlx>=0.22; platform_system == 'Darwin'", "pyopencl>=2025.2"]
 
 [project.scripts]
 fast-mlsirm = "fast_mlsirm.cli:main"

python/fast_mlsirm/backend.py

@@ -0,0 +1,49 @@
+from __future__ import annotations
+
+import importlib
+
+
+VALID_COMPUTE_BACKENDS = {"cpu", "cuda", "mlx", "opencl"}
+
+
+def normalize_backend(name: str) -> str:
+    backend = str(name).strip().lower()
+    if backend not in VALID_COMPUTE_BACKENDS:
+        raise ValueError(f"compute_backend must be one of {sorted(VALID_COMPUTE_BACKENDS)}")
+    return backend
+
+
+def ensure_backend_available(name: str) -> str:
+    backend = normalize_backend(name)
+    if backend == "cpu":
+        return backend
+    if backend == "cuda":
+        _require_module("cupy", "CUDA backend requires cupy.")
+        return backend
+    if backend == "mlx":
+        _require_module("mlx.core", "MLX backend requires mlx.")
+        return backend
+    _require_opencl()
+    return backend
+
+
+def _require_module(module_name: str, message: str) -> None:
+    try:
+        importlib.import_module(module_name)
+    except Exception as exc:  # pragma: no cover - depends on runtime environment
+        raise ValueError(message) from exc
+
+
+def _require_opencl() -> None:
+    try:
+        import pyopencl as cl
+    except Exception as exc:  # pragma: no cover - depends on runtime environment
+        raise ValueError("OpenCL backend requires pyopencl.") from exc
+
+    try:
+        platforms = cl.get_platforms()
+    except Exception as exc:  # pragma: no cover - depends on runtime environment
+        raise ValueError("OpenCL backend is unavailable because no OpenCL platform was found.") from exc
+
+    if not platforms:
+        raise ValueError("OpenCL backend is unavailable because no OpenCL platform was found.")

python/fast_mlsirm/cli.py

@@ -95,6 +95,12 @@ def main(argv: list[str] | None = None) -> int:
     fit_cmd.add_argument("--model", default="MLS2PLM", help="Model type to fit (default: MLS2PLM).")
     fit_cmd.add_argument("--latent-dim", type=int, default=2, help="Latent dimensionality for person traits (default: 2).")
     fit_cmd.add_argument("--optimizer", choices=["adam", "lbfgs", "adam_lbfgs"], default="adam_lbfgs", help="Optimizer to use (default: adam_lbfgs).")
+    fit_cmd.add_argument(
+        "--device",
+        choices=["cpu", "cuda", "mlx", "opencl"],
+        default="cpu",
+        help="Compute backend for fitting (default: cpu).",
+    )
     fit_cmd.add_argument("--max-iter", type=int, default=100, help="Maximum number of iterations for the optimizer (default: 100).")
     fit_cmd.add_argument("--n-restarts", type=int, default=1, help="Number of random restarts (default: 1).")
     fit_cmd.add_argument("--seed", type=int, default=1, help="Random seed for fitting (default: 1).")
@@ -126,6 +132,12 @@ def main(argv: list[str] | None = None) -> int:
     dim.add_argument("--folds", type=int, default=5, help="Number of validation folds (default: 5).")
     dim.add_argument("--model", default="MLS2PLM", help="Model type to fit (default: MLS2PLM).")
     dim.add_argument("--optimizer", choices=["adam", "lbfgs", "adam_lbfgs"], default="adam_lbfgs", help="Optimizer to use (default: adam_lbfgs).")
+    dim.add_argument(
+        "--device",
+        choices=["cpu", "cuda", "mlx", "opencl"],
+        default="cpu",
+        help="Compute backend for fitting (default: cpu).",
+    )
     dim.add_argument("--max-iter", type=int, default=100, help="Maximum iterations per fold fit (default: 100).")
     dim.add_argument("--n-restarts", type=int, default=1, help="Random restarts per fold fit (default: 1).")
     dim.add_argument("--seed", type=int, default=1, help="Random seed for folds and fitting (default: 1).")
@@ -280,6 +292,7 @@ def main(argv: list[str] | None = None) -> int:
                 max_iter=args.max_iter,
                 n_restarts=args.n_restarts,
                 seed=args.seed,
+                compute_backend=args.device,
             ),
         )
         save_dimensionality_diagnostics(diagnostics, args.out)
@@ -419,6 +432,7 @@ def main(argv: list[str] | None = None) -> int:
             max_iter=args.max_iter,
             n_restarts=args.n_restarts,
             seed=args.seed,
+            compute_backend=args.device,
         ),
     )
     save_fit_result(result, args.out)
@@ -431,6 +445,7 @@ def main(argv: list[str] | None = None) -> int:
             "out": str(args.out),
             "model": result.model,
             "optimizer": result.optimizer,
+            "compute_backend": result.compute_backend,
             "objective": float(result.objective),
             "convergence_status": result.convergence_status,
             "n_iter": int(result.n_iter),

python/fast_mlsirm/config.py

@@ -2,6 +2,7 @@
 
 from dataclasses import dataclass
 
+from .backend import normalize_backend
 
 VALID_MODELS = {"MIRT", "MLS2PLM", "MLSRM", "ULS2PLM", "ULSRM"}
 VALID_OPTIMIZERS = {"adam", "lbfgs", "adam_lbfgs"}
@@ -66,6 +67,7 @@ class FitConfig:
     gradient_clip: float | None = 100.0
     lbfgs_history: int = 10
     verbose: int = 0
+    compute_backend: str = "cpu"
     penalty: PenaltyConfig = PenaltyConfig()
 
     def normalized_model(self) -> str:
@@ -89,3 +91,4 @@ def validate(self) -> None:
             raise ValueError("init_gamma must be > 0")
         if self.eps_distance <= 0:
             raise ValueError("eps_distance must be > 0")
+        normalize_backend(self.compute_backend)

python/fast_mlsirm/fit.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 
+from .backend import ensure_backend_available
 from .config import FitConfig
 from .math import logit, normalize_latent_positions, standardize
 from .objective import model_flags, neg_loglik_and_grad, prepare_response, validate_factor_id
@@ -18,6 +19,7 @@ def fit(
 ) -> FitResult:
     config = config or FitConfig()
     config.validate()
+    compute_backend = ensure_backend_available(config.compute_backend)
     model = config.normalized_model()
 
     y, observed = prepare_response(responses, mask)
@@ -30,7 +32,7 @@ def fit(
 
     best: FitResult | None = None
     for restart in range(config.n_restarts):
-        candidate = _fit_single_restart(restart, config, y, observed, factors, n_dims, model)
+        candidate = _fit_single_restart(restart, config, y, observed, factors, n_dims, model, compute_backend)
         if best is None or candidate.objective < best.objective:
             best = candidate
 
@@ -47,11 +49,12 @@ def _fit_single_restart(
     factors: np.ndarray,
     n_dims: int,
     model: str,
+    compute_backend: str,
 ) -> FitResult:
     rng = np.random.default_rng(config.seed + restart)
     params0 = _initial_params(y, observed, factors, n_dims, config.latent_dim, config, rng)
     x0 = _pack(params0, model)
-    objective = _make_objective(y, observed, factors, params0, config)
+    objective = _make_objective(y, observed, factors, params0, config, compute_backend)
 
     x = x0
     obj_trace: list[float] = []
@@ -76,14 +79,17 @@ def _fit_single_restart(
     final_params = _unpack(x, params0, model)
     if model != "MIRT":
         final_params = normalize_latent_positions(final_params)
-    final_obj, _, final_loglik = neg_loglik_and_grad(y, factors, final_params, config, mask=observed)
+    final_obj, _, final_loglik = neg_loglik_and_grad(
+        y, factors, final_params, config, mask=observed, compute_backend=compute_backend
+    )
     obj_trace.append(final_obj)
     loglik_trace.append(final_loglik)
 
     candidate = FitResult(
         params=final_params,
         model=model,
         optimizer=config.optimizer,
+        compute_backend=compute_backend,
         objective=final_obj,
         loglik_trace=loglik_trace,
         objective_trace=obj_trace,
@@ -126,12 +132,13 @@ def _make_objective(
     factor_id: np.ndarray,
     template: MLSIRMParams,
     config: FitConfig,
+    compute_backend: str,
 ) -> Callable[[np.ndarray], tuple[float, np.ndarray, float]]:
     model = config.normalized_model()
 
     def objective(x: np.ndarray) -> tuple[float, np.ndarray, float]:
         params = _unpack(x, template, model)
-        obj, grad, loglik = neg_loglik_and_grad(y, factor_id, params, config, mask=observed)
+        obj, grad, loglik = neg_loglik_and_grad(y, factor_id, params, config, mask=observed, compute_backend=compute_backend)
         grad_vec = _pack(grad, model)
         if config.gradient_clip is not None:
             norm = float(np.linalg.norm(grad_vec))

python/fast_mlsirm/io.py

@@ -54,6 +54,7 @@ def save_fit_result(result: FitResult, run_dir: str | Path) -> None:
     summary = {
         "model": result.model,
         "optimizer": result.optimizer,
+        "compute_backend": result.compute_backend,
         "objective": result.objective,
         "convergence_status": result.convergence_status,
         "n_iter": result.n_iter,