"""Numerical BLIM trade-offs by parameter family (Stefanutti et al., 2012).
Independent linear-algebra implementation: singular-value decompositions,
not a port of the R/MATLAB row-reduction code. Null-space bases are not
unique, so compare their spans, ranks and prediction derivatives rather
than coefficients of a particular software's basis.
"""
from __future__ import annotations
from collections.abc import Mapping
from dataclasses import dataclass
from itertools import combinations
import numpy as np
from knowledgespaces._patterns import _positive_integer
from knowledgespaces.estimation.constraints import BLIMConstraints, _compile_constraints
from knowledgespaces.estimation.identifiability import _as_vector, blim_jacobian
from knowledgespaces.estimation.prediction import ItemParameter, StateParameter, _item_vector
from knowledgespaces.structures.knowledge_structure import KnowledgeStructure
[docs]
@dataclass(frozen=True)
class BLIMParameter:
"""One free Jacobian coordinate, including equality groups.
``family`` is beta, eta, or pi. ``items`` labels an error parameter's
equality class. For pi, ``state`` labels the coordinate; the mass of
the first canonical state is one minus the other masses. Fixed
parameters have no coordinate. Equality groups appear once.
"""
family: str
items: tuple[str, ...] = ()
state: frozenset[str] | None = None
[docs]
@dataclass(frozen=True)
class JacobianBlock:
"""Rank and orthonormal null-space basis for selected families.
``columns`` indexes the complete report's parameter list and Jacobian.
``null_space`` has one row per selected free coordinate and one column
per null direction. It includes dependencies within a family as well
as dependencies across families. Arrays are read-only. Multiplying
the selected Jacobian by this basis gives numerical zero, within the
reported absolute singular-value cutoff ``rank_tolerance``.
"""
families: tuple[str, ...]
columns: tuple[int, ...]
rank: int
rank_tolerance: float
singular_values: np.ndarray
null_space: np.ndarray
@property
def nullity(self) -> int:
"""Number of numerically null parameter directions."""
return len(self.columns) - self.rank
[docs]
@dataclass(frozen=True)
class BLIMTradeoffs:
"""Pointwise diagnostics for all seven beta/eta/pi block combinations.
``beta``, ``eta`` and ``pi`` record the actual evaluation point, after
constraints, in sorted item and canonical state order. ``jacobian``
uses lexicographic binary response order. Arrays are read-only.
Deficiency at one singular point is not proof of generic or global
non-identifiability. At boundaries, null directions describe the
polynomial derivative and need not be feasible probability directions.
Use check_identifiability for multiple-point rank screening.
"""
items: tuple[str, ...]
states: tuple[frozenset[str], ...]
beta: np.ndarray
eta: np.ndarray
pi: np.ndarray
parameters: tuple[BLIMParameter, ...]
jacobian: np.ndarray
blocks: tuple[JacobianBlock, ...]
[docs]
def block(self, *families: str) -> JacobianBlock:
"""Select a nonempty subset of beta, eta, pi, in any argument order."""
if len(families) == len(set(families)):
for block in self.blocks:
if set(block.families) == set(families):
return block
raise ValueError("Choose a nonempty subset of 'beta', 'eta', 'pi', without duplicates.")
[docs]
def blim_tradeoffs(
structure: KnowledgeStructure,
*,
beta: ItemParameter = 0.1,
eta: ItemParameter = 0.1,
pi: StateParameter = None,
constraints: BLIMConstraints | None = None,
rank_tolerance: float | None = None,
max_items: int | None = None,
max_memory_bytes: int = 512_000_000,
) -> BLIMTradeoffs:
"""Inspect ranks and null directions at one specified BLIM parameter point.
Uses the same constraints and coordinate convention as blim_jacobian.
``rank_tolerance`` is an optional nonnegative absolute SVD cutoff,
shared by all blocks. By default each block uses NumPy's convention
max(shape)*eps*largest_singular_value. Reported bases describe linear
dependencies, not unique parameter remedies or confidence intervals.
Every block is evaluated even if the full Jacobian has full rank.
The complete beta/eta/pi block handles dependencies involving three
families without assuming that pairwise ranks characterize them.
Extra memory for seven bases and SVD work is guarded by a conservative
allocation estimate; BLAS workspace and process RSS are not bounded.
"""
max_memory_bytes = _positive_integer(max_memory_bytes, "max_memory_bytes")
if rank_tolerance is not None and (not np.isfinite(rank_tolerance) or rank_tolerance < 0):
raise ValueError("rank_tolerance must be finite and nonnegative.")
items, states = sorted(structure.domain), list(structure)
compiled = _compile_constraints(constraints or BLIMConstraints(), items, states)
npar = compiled.n_parameters(len(states))
# Account for stored bases, Jacobian copies and SVD U/V, before enumeration.
extra_bytes = 8 * (10 * npar * npar + 4 * (2 ** len(items)) * npar)
if extra_bytes >= max_memory_bytes:
raise MemoryError("Trade-off diagnostics exceed max_memory_bytes.")
jacobian = blim_jacobian(
structure,
beta=beta,
eta=eta,
pi=pi,
constraints=constraints,
max_items=max_items,
max_memory_bytes=max_memory_bytes - extra_bytes,
)
resolver = _as_vector if constraints is None else _item_vector
b, e = resolver(beta, items, "beta"), resolver(eta, items, "eta")
compiled.beta.project(b, clip_free=False)
compiled.eta.project(e, clip_free=False)
if compiled.pi is not None:
prior = compiled.pi.copy()
elif pi is None:
prior = np.full(len(states), 1 / len(states))
elif isinstance(pi, Mapping):
prior = np.array([pi[state] for state in states], dtype=float)
else:
prior = np.array(pi, dtype=float, copy=True)
parameters = []
for family, spec in (("beta", compiled.beta), ("eta", compiled.eta)):
for group, fixed in zip(spec.groups, spec.fixed, strict=True):
if fixed is None:
parameters.append(BLIMParameter(family, tuple(items[i] for i in group)))
if compiled.pi is None:
parameters.extend(BLIMParameter("pi", state=state) for state in states[1:])
blocks = []
for size in range(1, 4):
for families in combinations(("beta", "eta", "pi"), size):
columns = tuple(i for i, p in enumerate(parameters) if p.family in families)
matrix = jacobian[:, columns]
if columns:
# Tall matrices need reduced U; wide ones need all right vectors.
_, singular, vh = np.linalg.svd(matrix, full_matrices=len(matrix) < len(columns))
cutoff = (
float(max(matrix.shape) * np.finfo(float).eps * singular[0])
if rank_tolerance is None
else rank_tolerance
)
rank = int(np.count_nonzero(singular > cutoff))
basis = vh[rank:].T.copy()
else:
singular, basis = np.empty(0), np.empty((0, 0))
rank, cutoff = 0, 0.0 if rank_tolerance is None else rank_tolerance
singular.flags.writeable = basis.flags.writeable = False
blocks.append(JacobianBlock(families, columns, rank, cutoff, singular, basis))
for array in (b, e, prior, jacobian):
array.flags.writeable = False
return BLIMTradeoffs(
tuple(items), tuple(states), b, e, prior, tuple(parameters), jacobian, tuple(blocks)
)