"""Random quasi orders and canonical BLIM simulations on their downsets.
The Bernoulli-pair/closure scheme is described by DAKS::simu (2.1-3).
It does not sample uniformly from quasi orders or control final density.
"""
from __future__ import annotations
from collections.abc import Collection, Mapping
from dataclasses import dataclass
import numpy as np
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation.blim_em import ResponseMatrix
from knowledgespaces.estimation.simulate import simulate_blim
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
from knowledgespaces.structures.relations import SurmiseRelation
from knowledgespaces.structures.set_family import SetFamily
def _domain(items: Collection[str], max_pairs: int) -> list[str]:
max_pairs = _positive_integer(max_pairs, "max_pairs")
labels = list(items)
if not labels or any(not isinstance(q, str) for q in labels) or len(set(labels)) != len(labels):
raise ValueError("items must contain distinct string labels on a nonempty domain.")
if len(labels) * (len(labels) - 1) > max_pairs:
raise MemoryError("Relation exceeds max_pairs; quadratic pair storage is required.")
return sorted(labels)
[docs]
def random_surmise_relation(
items: Collection[str],
delta: float,
*,
seed: int | np.random.Generator | None = None,
max_pairs: int = 1_000_000,
) -> SurmiseRelation:
"""Draw nonreflexive directed pairs independently, then close transitively.
Every ordered pair of distinct items has inclusion probability delta
before closure; reflexivity is implicit. Draws follow lexicographic
item order, row by row, omitting the diagonal. (a,b) means a is required
for b. Cycles produce equivalent items and are retained. delta must be
finite in [0,1]; 0 gives the identity, 1 the universal relation.
This is the DAKS documented sampling scheme, not a uniform law on
quasi orders. Transitive closure changes marginal edge probabilities.
max_pairs bounds the number of candidate nonreflexive pairs before
allocation; closure also has potentially substantial time cost.
"""
labels = _domain(items, max_pairs)
if not np.isfinite(delta) or not 0 <= delta <= 1:
raise ValueError("delta must be finite and in [0,1].")
rng = np.random.default_rng(seed)
pairs: list[tuple[str, str]] = []
for a in labels:
draws = rng.random(len(labels) - 1)
others = (b for b in labels if b != a)
pairs.extend((a, b) for b, u in zip(others, draws, strict=True) if u < delta)
return SurmiseRelation(labels, pairs).transitive_closure()
[docs]
@dataclass(frozen=True)
class QuasiOrdinalSimulation:
"""Generated relation, all compatible states, and respondent-aligned data.
data has one row per respondent, with no aggregation; true_states[r]
identifies respondent r's latent state. The structure is quasi-ordinal
but need not be a learning space when items are equivalent.
"""
data: ResponseMatrix
relation: SurmiseRelation
structure: KnowledgeStructure
true_states: tuple[frozenset[str], ...]
[docs]
def simulate_quasiordinal(
items: Collection[str],
n_respondents: int,
*,
relation: SurmiseRelation | None = None,
delta: float | None = None,
beta: float | Mapping[str, float] = 0.1,
eta: float | Mapping[str, float] = 0.1,
pi: Mapping[frozenset[str], float] | np.ndarray | None = None,
seed: int | np.random.Generator | None = None,
max_pairs: int = 1_000_000,
max_states: int = 100_000,
max_memory_bytes: int = 512_000_000,
) -> QuasiOrdinalSimulation:
"""Simulate a BLIM on a supplied or randomly generated quasi order.
Supply exactly one of relation (already transitive, matching items) or
delta (passed to random_surmise_relation). All downsets are enumerated
with a max_states bound. State masses are uniform by default; pi and
beta/eta follow simulate_blim, including heterogeneous error rates.
Output items are sorted; pi arrays follow canonical structure order.
Response probabilities are P(1|mastered)=1-beta and P(1|unmastered)=eta.
DAKS 2.1-3's sequential response mutation instead makes its effective
slip ce*(1-lg); that implementation discrepancy is not reproduced.
A shared local generator drives relation, latent-state and response
draws. The memory preflight estimates numerical working arrays; Python
set/object overhead is additionally controlled by max_pairs/max_states.
"""
labels = _domain(items, max_pairs)
n_respondents = _positive_integer(n_respondents, "n_respondents")
max_states = _positive_integer(max_states, "max_states")
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if (relation is None) == (delta is None):
raise ValueError("Supply exactly one of relation or delta.")
sample_bytes = 48 * n_respondents * len(labels)
if sample_bytes > max_memory_bytes:
raise MemoryError("Simulation response arrays exceed max_memory_bytes.")
rng = np.random.default_rng(seed)
if relation is None:
assert delta is not None
relation = random_surmise_relation(labels, delta, seed=rng, max_pairs=max_pairs)
elif relation.items != set(labels):
raise ValueError("Relation domain must match items exactly.")
elif relation.transitive_closure() != relation:
raise ValueError("Supply a transitively closed relation, not its generators.")
principal = [{q} | relation.prerequisites_of(q) for q in labels]
family = SetFamily(labels, principal).union_closure(max_sets=max_states)
if len(family) > max_states:
raise ValueError("Structure exceeds max_states.")
if sample_bytes + 16 * len(family) * len(labels) > max_memory_bytes:
raise MemoryError("Simulation state arrays exceed max_memory_bytes.")
structure = family.to_knowledge_structure()
sample = simulate_blim(
structure, n_respondents, beta=beta, eta=eta, pi=pi, seed=rng, return_states=True
)
assert isinstance(sample, tuple)
data, states = sample
return QuasiOrdinalSimulation(data, relation, structure, tuple(states))