"""Reproducible random-start search for observed-data BLIM likelihoods."""
from __future__ import annotations
import warnings
from dataclasses import dataclass, replace
from typing import Literal
import numpy as np
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation._starts import (
_random_error_starts,
_random_prior_start,
_validate_start_options,
)
from knowledgespaces.estimation.blim_em import ConvergenceWarning, _warn_degenerate
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.incomplete import (
IncompleteBLIMEstimate,
IncompleteResponseMatrix,
UnobservedItemWarning,
_validate_incomplete_fit,
estimate_blim_incomplete,
)
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
[docs]
@dataclass(frozen=True)
class IncompleteBLIMRestart:
"""One attempted fit, including failures and runs reaching the iteration cap.
Effective starts include fixed/equality projection, in the returned fit's
item and canonical state orders. ``objective`` is the observed-data log
likelihood, equal to ``log_likelihood``. An errored run has NaN scores,
zero completed iterations and a nonempty ``error``. Nonconvergence is
separate: a capped fit has finite scores and is eligible for selection.
"""
beta_init: tuple[float, ...]
eta_init: tuple[float, ...]
pi_init: tuple[float, ...]
log_likelihood: float
objective: float
n_iterations: int
converged: bool
error: str | None = None
[docs]
class IncompleteBLIMRestartError(RuntimeError):
"""Every start failed; attempted starts and errors remain in ``restarts``."""
def __init__(self, restarts: tuple[IncompleteBLIMRestart, ...]) -> None:
self.restarts = restarts
super().__init__(f"All {len(restarts)} incomplete BLIM starts failed.")
[docs]
def estimate_blim_incomplete_restarts(
structure: KnowledgeStructure,
data: IncompleteResponseMatrix,
*,
n_restarts: int = 10,
seed: int | np.random.Generator | None = None,
init_range: tuple[float, float] = (0.01, 0.4),
init_strategy: Literal["uniform", "pks"] = "uniform",
constraints: BLIMConstraints | None = None,
max_iter: int = 5000,
tol: float = 1e-7,
max_memory_bytes: int = 512_000_000,
) -> IncompleteBLIMEstimate:
"""Maximize observed-data BLIM likelihood across random initializations.
``uniform`` draws independent error parameters in ``init_range``, halves
paired values until beta+eta<.95, and starts pi at equal masses. ``pks``
draws errors from U(0,1) with sequential beta/eta reflection and pi from
Dirichlet(1, ..., 1), the uniform-simplex law of pks 0.7-0. Constraints
override starts; the fit does not impose beta+eta<1 on subsequent iterates.
The name describes initialization only, not an incomplete-data pks fit.
Every run retains effective starts, likelihood, objective, convergence,
iterations and numerical failure. Select the largest finite likelihood;
first run wins exact ties, even if unfinished. A selected unfinished fit
warns; all-error searches raise IncompleteBLIMRestartError with records.
Invalid arguments and resource limits raise directly, not as fit failures.
Seed/Generator controls a local stream, never NumPy's global RNG.
This changes optimization, not missingness assumptions, identifiability,
or the absence of a global-optimum guarantee. See estimate_blim_incomplete.
"""
n_restarts = _positive_integer(n_restarts, "n_restarts")
_validate_start_options(init_range, init_strategy)
max_iter, max_memory_bytes = _validate_incomplete_fit(
structure, data, max_iter, tol, max_memory_bytes
)
compiled = _compile_constraints(
constraints or BLIMConstraints(), list(data.items), list(structure)
)
if 8 * n_restarts * (2 * data.n_items + len(structure) + 6) > max_memory_bytes:
raise MemoryError("Incomplete BLIM restart records exceed max_memory_bytes.")
rng = np.random.default_rng(seed)
records: list[IncompleteBLIMRestart] = []
best: IncompleteBLIMEstimate | None = None
selected = 0
for index in range(n_restarts):
beta, eta = _random_error_starts(rng, data.n_items, init_range, init_strategy)
prior = _random_prior_start(rng, len(structure), init_strategy)
compiled.beta.project(beta)
compiled.eta.project(eta)
if compiled.pi is not None:
prior = compiled.pi.copy()
error = None
fit = None
try:
with warnings.catch_warnings():
warnings.simplefilter("ignore", ConvergenceWarning)
if index:
warnings.simplefilter("ignore", UnobservedItemWarning)
fit = estimate_blim_incomplete(
structure,
data,
beta_init=beta,
eta_init=eta,
pi_init=prior,
constraints=constraints,
max_iter=max_iter,
tol=tol,
max_memory_bytes=max_memory_bytes,
)
except ValueError as exc:
if str(exc) != "Observed responses are impossible under the fixed/initial parameters.":
raise
error = str(exc)
except FloatingPointError as exc:
error = str(exc)
if fit is not None and not np.isfinite(fit.log_likelihood):
error = "Non-finite observed-data log likelihood."
fit = None
records.append(
IncompleteBLIMRestart(
tuple(beta.tolist()),
tuple(eta.tolist()),
tuple(prior.tolist()),
float("nan") if fit is None else fit.log_likelihood,
float("nan") if fit is None else fit.objective_history[-1],
0 if fit is None else fit.n_iterations,
False if fit is None else fit.converged,
error,
)
)
if fit is not None and (best is None or fit.log_likelihood > best.log_likelihood):
best, selected = fit, index
if best is None:
raise IncompleteBLIMRestartError(tuple(records))
if not best.converged or any(r.error is not None for r in records):
failed = sum(r.error is not None for r in records)
capped = sum(r.error is None and not r.converged for r in records)
warnings.warn(
f"Incomplete BLIM search: {failed} failed, {capped} capped starts; "
f"selected start {selected} converged={best.converged}.",
ConvergenceWarning,
stacklevel=2,
)
_warn_degenerate(best.degenerate_items)
return replace(best, restarts=tuple(records), selected_restart=selected)