"""Scientific BLIM/SLM summaries with explicit data and sample conventions."""
from __future__ import annotations
import math
from dataclasses import asdict, dataclass
import numpy as np
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation._data_signature import _data_signature
from knowledgespaces.estimation.blim_em import BLIMEstimate, GoodnessOfFit, ResponseMatrix
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.incomplete import IncompleteBLIMEstimate
from knowledgespaces.estimation.prediction import (
_constraint_feasibility,
_item_vector,
_state_matrix,
_state_vector,
)
from knowledgespaces.estimation.slm import SLMEstimate
from knowledgespaces.estimation.tradeoffs import BLIMTradeoffs, blim_tradeoffs
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
Fit = BLIMEstimate | IncompleteBLIMEstimate
[docs]
@dataclass(frozen=True)
class MinimumDiscrepancy:
"""Hamming minima over states feasible under declared fixed-zero errors.
``distances`` retains input row order; infeasible rows have +inf.
``counts[d]`` is the supplied weight at finite distance d (0..|Q|).
``infeasible_weight`` is separate, so no observations silently vanish.
``mean`` is +inf if any positive-weight row is infeasible. This is a
distance distribution, not a BLIM likelihood or state-posterior summary.
``training_data`` compares labeled empirical measures, not participant IDs;
None means the fit has no training fingerprint. Arrays are read-only.
"""
items: tuple[str, ...]
distances: np.ndarray
counts: np.ndarray
total_weight: float
infeasible_weight: float
mean: float
training_data: bool | None
[docs]
@dataclass(frozen=True)
class Coefficient:
"""One labeled reported coefficient; equality groups appear only once.
Pi rows contain actual state masses (not free simplex coordinates).
In SLM they are ``derived=True`` from g, rather than freely estimated.
Fixed error groups use the merged restrictions, including propagation
of a fixed value through overlapping equality groups.
"""
family: str
value: float
items: tuple[str, ...] = ()
state: frozenset[str] | None = None
fixed: bool = False
derived: bool = False
[docs]
@dataclass(frozen=True)
class BLIMReport:
"""Immutable numerical report with separate training and supplied-data roles.
Mastery is sum_K pi[K] 1(q in K). Expected slips and guesses are,
respectively, beta[q]*mastery[q] and eta[q]*(1-mastery[q]): marginal
probabilities of error for a new respondent taking the entire test.
Their sums are expected error counts per complete test, not fitted
errors conditional on observed responses or missingness masks.
AIC/BIC and ``training_gof`` belong to the stored fit, irrespective of
the data supplied for discrepancy. The nominal parameter count is not
an identified dimension. AICc is a conventional heuristic, with no
general validity guarantee for singular/boundary BLIMs or for MD fits.
For missing data, no automatic GOF or discrepancy is supplied.
"""
model: str
method: str
items: tuple[str, ...]
states: tuple[frozenset[str], ...]
converged: bool
n_iterations: int
n_parameters: int
coefficients: tuple[Coefficient, ...]
constraints: BLIMConstraints
mastery: np.ndarray
expected_slips: np.ndarray
expected_guesses: np.ndarray
training_log_likelihood: float
training_gof: GoodnessOfFit | None
AIC: float
BIC: float
BIC_npatterns: float | None
AICc: float | None
aicc_sample_size: float | None
aicc_sample_convention: str | None
aicc_unavailable_reason: str | None
discrepancy: MinimumDiscrepancy | None
tradeoffs: BLIMTradeoffs | None
@property
def expected_total_errors(self) -> float:
"""Expected number of slips plus guesses per new complete test."""
return float(self.expected_slips.sum() + self.expected_guesses.sum())
[docs]
def to_dict(self) -> dict:
"""Return a detached JSON-ready export, with labels and interpretation.
Nonfinite numbers become None (JSON null), with infeasible weight
and AICc's unavailable reason retained. ``state`` is a sorted label
list, never a lossy string representation. No file is written.
"""
restrictions = {
"beta_fixed": dict(self.constraints.beta_fixed),
"eta_fixed": dict(self.constraints.eta_fixed),
"beta_equal": [list(g) for g in self.constraints.beta_equal],
"eta_equal": [list(g) for g in self.constraints.eta_equal],
"pi_fixed": None
if self.constraints.pi_fixed is None
else [
{"state": sorted(state), "value": value}
for state, value in self.constraints.pi_fixed.items()
],
}
result: dict = {
"model": self.model,
"method": self.method,
"items": list(self.items),
"states": [sorted(s) for s in self.states],
"converged": self.converged,
"n_iterations": self.n_iterations,
"n_parameters": self.n_parameters,
"coefficients": [
{
"family": c.family,
"value": c.value,
"items": list(c.items),
"state": None if c.state is None else sorted(c.state),
"fixed": c.fixed,
"derived": c.derived,
}
for c in self.coefficients
],
"constraints": restrictions,
"population_full_test": {
"mastery": dict(zip(self.items, self.mastery.tolist(), strict=True)),
"expected_slips": dict(zip(self.items, self.expected_slips.tolist(), strict=True)),
"expected_guesses": dict(
zip(self.items, self.expected_guesses.tolist(), strict=True)
),
"expected_total_errors": self.expected_total_errors,
},
"training": {
"log_likelihood": self.training_log_likelihood,
"gof": None if self.training_gof is None else asdict(self.training_gof),
"AIC": self.AIC,
"BIC": self.BIC,
"BIC_npatterns": self.BIC_npatterns,
"AICc": self.AICc,
"aicc_sample_size": self.aicc_sample_size,
"aicc_sample_convention": self.aicc_sample_convention,
"aicc_unavailable_reason": self.aicc_unavailable_reason,
},
"supplied_data_discrepancy": None
if self.discrepancy is None
else {
"items": list(self.discrepancy.items),
"distances": self.discrepancy.distances.tolist(),
"counts": self.discrepancy.counts.tolist(),
"distance_values": list(range(len(self.discrepancy.counts))),
"total_weight": self.discrepancy.total_weight,
"infeasible_weight": self.discrepancy.infeasible_weight,
"mean": self.discrepancy.mean,
"training_data": self.discrepancy.training_data,
},
"tradeoffs": None,
}
if self.tradeoffs is not None:
report = self.tradeoffs
result["tradeoffs"] = {
"scope": "Full-response BLIM derivative at the fitted parameters; not global identifiability or a missing-mask derivative.",
"items": list(report.items),
"states": [sorted(s) for s in report.states],
"parameters": [
{
"family": p.family,
"items": list(p.items),
"state": None if p.state is None else sorted(p.state),
}
for p in report.parameters
],
"jacobian": report.jacobian.tolist(),
"blocks": [
{
"families": list(b.families),
"columns": list(b.columns),
"rank": b.rank,
"nullity": b.nullity,
"rank_tolerance": b.rank_tolerance,
"singular_values": b.singular_values.tolist(),
"null_space": b.null_space.tolist(),
}
for b in report.blocks
],
}
return _finite_json(result)
def _finite_json(value):
if isinstance(value, dict):
return {key: _finite_json(item) for key, item in value.items()}
if isinstance(value, list):
return [_finite_json(item) for item in value]
if isinstance(value, float) and not math.isfinite(value):
return None
return value
[docs]
def aicc(log_likelihood: float, n_parameters: int, n_observations: float) -> float:
"""Conventional AIC + 2*k*(k+1)/(N-k-1), requiring an explicit N.
Returns NaN when N <= k+1: the correction is undefined there, rather
than a negative reward. N must be finite and positive, k a nonnegative
integer, and log likelihood finite. The MATLAB KST-toolbox uses total
respondent count for N and *negative* log likelihood internally. This
function takes the ordinary log likelihood. It is an algebraic criterion,
not a general small-sample bias correction for singular BLIM models.
Fractional N is accepted as an explicitly chosen descriptive convention;
it does not turn survey/analysis weights into independent observations.
"""
if isinstance(n_parameters, bool) or not isinstance(n_parameters, int) or n_parameters < 0:
raise ValueError("n_parameters must be a nonnegative integer.")
if isinstance(n_observations, bool) or not np.isfinite(n_observations) or n_observations <= 0:
raise ValueError("n_observations must be finite and positive.")
if not np.isfinite(log_likelihood):
raise ValueError("log_likelihood must be finite.")
if n_observations <= n_parameters + 1:
return math.nan
return float(
-2 * log_likelihood
+ 2 * n_parameters
+ 2 * n_parameters * (n_parameters + 1) / (n_observations - n_parameters - 1)
)
[docs]
def minimum_discrepancy(
estimate: Fit,
data: ResponseMatrix,
*,
chunk_size: int = 1024,
max_memory_bytes: int = 512_000_000,
) -> MinimumDiscrepancy:
"""Weighted minimum Hamming distances on explicitly complete responses.
Matches the defining distance in pks getMD: declared fixed-zero slip or
guess groups exclude impossible state/response pairs before minimizing.
An estimated numerical zero or zero prior mass does not exclude a state.
The excess inclusion radius and hyperbolic MD weights do not change the
minimum itself. No powerset of responses is enumerated; inputs may be
compressed frequencies or uncompressed observations, including zero rows.
"""
if not isinstance(estimate, (BLIMEstimate, IncompleteBLIMEstimate)):
raise TypeError("estimate must be a BLIM, SLM or incomplete BLIM estimate.")
if not isinstance(data, ResponseMatrix):
raise TypeError("minimum_discrepancy requires complete ResponseMatrix data.")
# ResponseMatrix is public and its arrays/lists are mutable. Revalidate
# a detached snapshot before distances or weights enter the report.
data = ResponseMatrix(
list(data.items),
np.array(data.patterns, copy=True),
None if data.counts is None else np.array(data.counts, copy=True),
)
if set(data.items) != set(estimate.items):
raise ValueError("data.items must match the estimate's domain.")
chunk_size = _positive_integer(chunk_size, "chunk_size")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
items, states = list(data.items), list(estimate.states)
b = _item_vector(estimate.beta_dict(), items, "beta")
e = _item_vector(estimate.eta_dict(), items, "eta")
prior = _state_vector(estimate.pi_dict(), states)
n, m, q = data.n_patterns, len(states), len(items)
if 8 * (7 * min(n, chunk_size) * m + 3 * m * q + n + q) > max_memory_bytes:
raise MemoryError("Minimum-discrepancy report exceeds max_memory_bytes.")
membership = _state_matrix(items, states)
distances = np.empty(n)
for start in range(0, n, chunk_size):
patterns = data.patterns[start : start + chunk_size].astype(float)
matrix = patterns @ (1 - membership).T + (1 - patterns) @ membership.T
feasible = _constraint_feasibility(
estimate.constraints, items, states, patterns, membership, b, e, prior
)
if feasible is not None:
matrix[~feasible] = np.inf
distances[start : start + len(patterns)] = matrix.min(axis=1)
finite = np.isfinite(distances)
counts = np.bincount(
distances[finite].astype(int), weights=data.effective_counts[finite], minlength=q + 1
)
infeasible_weight = float(data.effective_counts[~finite].sum())
mean = (
math.inf if infeasible_weight > 0 else float(counts @ np.arange(q + 1) / data.n_respondents)
)
signature = getattr(estimate, "data_signature", None)
training = None if signature is None else signature == _data_signature(data)
distances.flags.writeable = counts.flags.writeable = False
return MinimumDiscrepancy(
tuple(items), distances, counts, data.n_respondents, infeasible_weight, mean, training
)
[docs]
def blim_report(
estimate: Fit,
data: ResponseMatrix | None = None,
*,
aicc_sample_size: float | None = None,
include_tradeoffs: bool = False,
max_memory_bytes: int = 512_000_000,
) -> BLIMReport:
"""Summarize fitted BLIM/SLM coefficients, errors, mastery and diagnostics.
Optional complete ``data`` supplies only the distance distribution.
Stored GOF/criteria remain training quantities. AICc automatically uses
N only for fingerprint-matched training data with integer frequencies.
For an incomplete fit, the report cannot verify whether stored total
weights represent respondent counts, so supply ``aicc_sample_size``
explicitly. Held-out N is never substituted.
The explicit override is a user-selected *training* N convention.
Trade-offs are opt-in because they enumerate all complete patterns.
They describe the full-response BLIM derivative, also when the estimate
used missing data. For SLM use ``estimate.jacobian()`` instead: its state
masses are derived from g, so requesting BLIM trade-offs raises.
"""
if not isinstance(estimate, (BLIMEstimate, IncompleteBLIMEstimate)):
raise TypeError("estimate must be a BLIM, SLM or incomplete BLIM estimate.")
if not isinstance(include_tradeoffs, bool):
raise ValueError("include_tradeoffs must be a bool.")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if include_tradeoffs and isinstance(estimate, SLMEstimate):
raise ValueError(
"SLM trade-offs require its own parameterization; use estimate.jacobian()."
)
items, states = list(estimate.items), list(estimate.states)
b = _item_vector(estimate.beta, items, "beta")
e = _item_vector(estimate.eta, items, "eta")
prior = _state_vector(estimate.pi, states)
constraints = estimate.constraints or BLIMConstraints()
compiled = _compile_constraints(constraints, items, states)
if 8 * (3 * len(items) * len(states) + 8 * len(items) + 4 * len(states)) > max_memory_bytes:
raise MemoryError("Statistical report exceeds max_memory_bytes.")
membership = _state_matrix(items, states)
_constraint_feasibility(
constraints, items, states, np.empty((0, len(items))), membership, b, e, prior
)
mastery = prior @ membership
slips, guesses = b * mastery, e * (1 - mastery)
coefficients = []
for name, values, spec in (("beta", b, compiled.beta), ("eta", e, compiled.eta)):
for group, fixed in zip(spec.groups, spec.fixed, strict=True):
coefficients.append(
Coefficient(
name,
float(values[group[0]]),
tuple(items[i] for i in group),
fixed=fixed is not None,
)
)
coefficients.extend(
Coefficient(
"pi",
float(p),
state=s,
fixed=compiled.pi is not None,
derived=isinstance(estimate, SLMEstimate),
)
for s, p in zip(states, prior, strict=True)
)
if isinstance(estimate, SLMEstimate):
coefficients.extend(
Coefficient("g", float(g), (q,)) for q, g in zip(items, estimate.g, strict=True)
)
discrepancy = (
None
if data is None
else minimum_discrepancy(estimate, data, max_memory_bytes=max_memory_bytes)
)
incomplete = isinstance(estimate, IncompleteBLIMEstimate)
npar = estimate.npar if isinstance(estimate, IncompleteBLIMEstimate) else estimate.gof.npar
AIC, BIC = (
(estimate.AIC, estimate.BIC)
if isinstance(estimate, IncompleteBLIMEstimate)
else (estimate.gof.AIC, estimate.gof.BIC)
)
n = aicc_sample_size
convention = "explicit_user_training_sample_size" if n is not None else None
reason = None
if n is None and data is not None and discrepancy is not None and discrepancy.training_data:
if np.equal(data.effective_counts, np.floor(data.effective_counts)).all():
n, convention = data.n_respondents, "verified_training_respondent_count"
else:
reason = "Fractional weights require an explicit sample-size convention."
if n is None:
correction = None
reason = (
reason
or "Training sample size unavailable; supply it explicitly or provide verified complete training data."
)
else:
value = aicc(estimate.log_likelihood, npar, n)
correction = value if np.isfinite(value) else None
if correction is None:
reason = "AICc is undefined because N <= k + 1."
tradeoffs = None
if include_tradeoffs:
tradeoffs = blim_tradeoffs(
KnowledgeStructure(items, states),
beta=estimate.beta_dict(),
eta=estimate.eta_dict(),
pi=estimate.pi_dict(),
constraints=constraints,
max_memory_bytes=max_memory_bytes,
)
for array in (mastery, slips, guesses):
array.flags.writeable = False
return BLIMReport(
"incomplete BLIM" if incomplete else "SLM" if isinstance(estimate, SLMEstimate) else "BLIM",
"ML" if isinstance(estimate, IncompleteBLIMEstimate) else estimate.method,
tuple(items),
tuple(states),
estimate.converged,
estimate.n_iterations,
npar,
tuple(coefficients),
constraints,
mastery,
slips,
guesses,
estimate.log_likelihood,
None if isinstance(estimate, IncompleteBLIMEstimate) else estimate.gof,
AIC,
BIC,
None if isinstance(estimate, IncompleteBLIMEstimate) else estimate.gof.BIC_npatterns,
correction,
n,
convention,
reason,
discrepancy,
tradeoffs,
)